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Abstract 

We describe the input data and installation procedures of the code hfodd (v2.40h). The 
present write-up contains complete and comprehensive information that has originally been 
given in six independent publications. It is enhanced by the subject index and indexes of 
variables, input-data keywords, subroutines, and files that are used in this user guide. 
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1. INTRODUCTION 

The code hfodd (v2.40h) solves the nuclear Skyrme-Hartree-Fock or Skyrme-Hartree-Fock- 
Bogolyubov problem by using the Cartesian deformed harmonic-oscillator basis. It is a re- 
sult of a long-term project that is developed at the Institute of Theoretical Physics, Univer- 
sity of Warsaw and in many other collaborating institutions. The code is one of the main 



solvers that are developed within the UNEDF collaboration (http://unedf.org/). The code 



has been published and is available from the Computer Physics Communications Program Li- 



brary (http://www.cpc.cs.qub.ac.uk/cpc/). The description of the code, and how to use it, 
has been published in the Computer Physics Communications in the form of six independent 
publications, which are below referred to as I [1], II [2], III [3], IV [1], V [5], and VI [6]. 

The present user's guide summarizes articles I- VI and provides practical details presented 
in the form that is suitable for a rapid reference on how to run the code. It also contains details 
updated to those corresponding to version (v2.40h). Information pertaining to the physics 
issues and numerical methods are not repeated here, and the user is referred to the original 
publications and to the references given therein. 

Together with consulting the present user's guide, the reader is invited to visit the hfodd 



home page (http://www.fuw.edu.pl/~dobaczew/hfodd/hfodd.html), where many useful files 



are available for downloading and which contains information on future versions and extensions 
being developed. 

The present user's guide overrides the previous version published in Ref. [7j, which is now 
obsolete. 
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2. INPUT DATA FILE 

The code reads the input data from the standard FORTRAN input file, and a general structure 
of the input data file is defined by the following rules: 

1. Input data file is an ASCII file composed of independent input items. 

2. An item is composed of the keyword line, and of the data line which follows immediately 
the keyword hne. Two items contain no data lines, see point [6] below. 

3. The keyword hne contains in columns 1 through 10 a keyword which is a specific text 
defining the item. If the keyword has less than 10 characters it has to be obligatorily 
padded with appropriate number of trailing spaces. Columns beyond 10 are ignored and 
can be used to place arbitrary comments or texts. 

4. The data hne contains a prescribed number of data values which are read in the free 
FORTRAN FORMAT. The type of data (REAL, INTEGER, or CHARACTER) should match 
the definition of the given data line. The CHARACTER data should be provided starting 
from the 13-th column of the data line. 

5. Items can be separated by an arbitrary numbers of comment lines. A line is treated as 
a comment and ignored provided it does not contain in columns 1 through 10 any of the 
valid keywords. 

6. Two items contain no data lines. The first one is defined by the keyword EXECUTE, and 
starts a calculation for the currently defined input parameters. The second one is defined 
by the keyword ALL_DONE, and terminates the program. The latter item is not required 
if the compiler is able to properly recognize the end of input data file. When the end of 
input data file is encountered, the program terminates as if the keyword ALL_DONE was 
found. 

7. A given item may occur any number of times in different places of the input data file. 
Only the last one of the same items occurring before a given EXECUTE keyword is taken 
into account. 

8. The items appearing between two consecutive items EXECUTE define the data set, i.e, the 
set of input parameters, for which the calculation will start at the moment when the 
latter of the two such items is found. Between two consecutive items EXECUTE, the order 
of other items is arbitrary. 

9. All input parameters have the default values defined in the code (subroutine predef). 
Therefore, if the input data file contains only the item EXECUTE, and no other lines, a 
calculation for the default values is performed. In the present version, this results in 
performing 50 iterations for the superdeformed state of ^^^Dy at the angular frequency of 
hiu=0.5MeV. 

10. Only those items which define the values of input parameters which are different from the 
default values have to be included in the input data file. The values of input parameters 
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defined by any item stay in effect till another occurrence of the same item, or till the code 
terminates. This rule is valid irrespective of how many EXECUTE items follow the given 
item. 

In the following subsections we list all the valid keywords and their corresponding data lines. 
Titles of the subsections are introduced only to organize the keywords into groups pertaining 
to specific categories, but are otherwise irrelevant. The default values are given in data lines 
following the keyword lines. In the data lines we also show names of the corresponding variables, 
which are referred to in the text. 
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2. INPUT DATA FILE 



2.1. General data 

This section lists keywords pertaining to definitions of particle numbers, iterations, and phase 
space. 



2.1.1. Keyword: NUCLIDE 

86, 66 = IN_FIX, IZ_FIX 

Numbers of neutrons (IN_FIX) and protons (IZ_FIX) in the nucleus under consideration. Calcu- 
lations for odd and odd-odd nuclei can be done only either with IR0TAT=1 or with IPAIRI=1. 
In the latter case, the odd and odd-odd nuclei are described within the false-vacuum approxi- 
mation unless the quasiparticle blocking approximation is used for IPAHFB=1, see Sees. I2.7H2.7I 
For IPABCS=1, the blocking approximation is not yet implemented. For IN_FIX=0 or IZ_FIX=0, 
calculations are performed for one type of particles only, which allows to treat the atomic gases. 



2.1.2. Keyword: ITERATIONS 

50 = NO ITER 

For NOITER>0, the specified number of iterations is performed. Unless IC0NTI=1 or IREAWS=1, 
the code starts with iteration number by constructing the initial mean fields based on the 
Nilsson potential. Specifying NOITER=0 requests only this initial phase. If the iteration is 
restarted from the previously recorded potentials (IC0NTI=1), then the counting of iterations 
continues from the number of iterations performed in the previous run. Along with requesting 
the initial Djj^ transformations (INIINV>0 or INIKAR>0) or initial rotation (INIR0T=1) one 
must set N0ITER=1. 



2.1.3. Keyword: ITERAT_EPS 

0.0 = EPSITE 



2.1.4. Keyword: MAXANTIOSC 

1 = NULAST 

These two parameters govern the termination of the HF iterations according to the achieved 
stability of solutions. The stability of the HF energy has been defined in (1-37) as the difference 
between the total energies calculated from the single-particle energies and from the Skyrme 
functional. The HF iterations continue until the absolute value of the stability is smaller than 
EPSITE (in MeV) over NULAST consecutive iterations. When this condition is fulfilled, iteration 
procedure terminates and the final results are printed. This allows for an automated adjustment 
of the number of iterations that are required to achieve a given level of convergence. The number 
of iterations NOITER, can now be set to a large value, at which the iterations terminate if a 
stable solution is not found. 
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The default value of EPSITE=0.0 deactivates this option, and then the code continues up 
to NOITER iterations regardless of the achieved stability. If a non-zero value of EPSITE is used, 
a non-zero value of NULAST has to be used too. In practice, a value of NULAST=5 prevents 
an accidental termination of iterations in all cases when the stability energy gradually goes 
through zero, but the solution is not yet self-consistent. 



2.1.5. Keyword: PING-PONG 



0.0, 3 = EPSPNG, NUPING 



The code is able to detect the "ping-pong" divergence described in Sec. III-2.6, i.e., the situation 
when the HF iteration procedure gives oscillating results in every second iteration. Upon 
continuing the iteration, both sequences of results, i.e., those which correspond to the iteration 
numbers being even and odd, stay different but perfectly stable, and hence the correct self- 
consistent solution is never attained. 

The code recognizes such a situation by calculating the averages and variances of the stability 
energy (1-37), separately in the even and in the odd sequences of results, over the last NUPING 
pairs of iterations. The "ping-pong" divergence condition occurs when both variances become 
a factor EPSPNG smaller than the absolute value of the difference of the corresponding averages, 
i.e., when 



where 



A(5£)evcn < EPSPNG X l^^evcn " 5^odd|, 
A{5S)odd < EPSPNG X |5^even -5^odd|, 



/NUPING, 





/ NUPING- 


-1 
















/ NUPING- 


-1 


<^^odd = j 




6£n- 




V P=0 





2p-l 



/NUPING, 



(la) 
(lb) 



(2a) 
(2b) 



and 



even I 



A{6S) 



odd 



(NUPING-1 
{68n-2p - SS 
p=0 

(NUPING-1 
p=0 



1/2 



/NUPING, 

1/2 

/NUPING. 



(3a) 



(3b) 



Here, n denotes the number of the last accomplished HF iteration, and denotes the stability 
energy (1-37) obtained in the n'-th iteration. 

The default value of EPSPNG=0.0 deactivates this option, and then the code continues up 
to NOITER iterations regardless of the "ping-pong" divergence. If a non-zero value of EPSPNG 
is used, a value of NUPING>1 has to be used too. In practice, values of EPSPNG=0.01 and 
NUPING=3 allow for an efficient detection of the "ping-pong" divergence condition. 
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2. INPUT DATA FILE 



Upon discovering the "ping-pong" divergence, the HF iterations are terminated and a table 
of absolute values of maximum differences of single-particle observables between the two se- 
quences of iterations is printed, see Sec. III-2.6. These maximum differences are determined for 
states in each of the charge-parity-signature, charge-simplex, parity, or charge blocks, depend- 
ing on the conserved symmetries, and separately for particle and hole states. Whenever such 
a maximum difference is found for a particle state and for a hole state with adjacent indices, 
such a pair is proposed as a candidate for the diabatic blocking calculation. 



2.1.6. Keyword: CHAOTIC 

= NUCHAO 

The code is able to detect the "chaotic" divergence which occurs when the HF iterations give 
results which chaotically vary from one iteration to another one. The code recognizes such 
a divergence by finding the local maxima Mfc, k=l,2,. . . , in the sequence of absolute values 
of the stability energies (1-37), obtained in the entire series of the HF iterations performed. 
The "chaotic" divergence condition occurs when the code finds NUCHAO positive differences 
Mfc— Mfc_i. When this condition occurs, iteration procedure terminates and the final results 
are printed. 

For NUCHAO=0 (the default value) the code does not check whether the "chaotic" divergence 
occurs or not. In practice, a value of NUCHA0=5 allows for an efficient detection of the "chaotic" 
divergence condition. However, for a small value of NUCHAO and a small value of EPSPNG, the 
"ping-pong" divergence can sometimes be mistaken for the "chaotic" divergence. If one is inter- 
ested in the diabatic-blocking data, printed after the "ping-pong" divergence, the recommended 
value of NUCHA0=5 should be increased to 10 or more. 



2.1.7. Keyword: PHASESPACE 

0, 0, 0, =NUMBSP(0,0), NUMBSP(1,0), 
NUMBSP(0,1), NUMBSP(1,1) 

Numbers of the lowest mean-field eigenstates which are kept after the diagonalization of the 
mean-field Hamiltonians in the four charge-simplex blocks: {s,q) = {+i,n), {—i,n), {+i,p), 
{—i,p). All other eigenstates are discarded. If any of these numbers is equal to zero (the 
default value), the code sets it equal to the number of neutrons IN_FIX (for q=n) or protons 
IZ_FIX (for q=p). When no symmetry is conserved, the numbers of states kept in charge blocks 
are determined by sums NUMBSP(0,0) + NUMBSP(1,0) and NUMBSP(0,1) + NUMBSP(1,1) 

For calculations without pairing, i.e. for IPAIRI=0, the user is responsible for using values of 
NUMBSP large enough to accommodate all wave functions which might be useful for the required 
vacuum and particle-hole configurations, see Sec, II-3.4. In practice, the use of the default 
values described above is recommended as a safe option. The size of the matrices defined by 
the NDSTAT parameter can be reduced by the user when small values of NUMBSP are used. 

For calculations with the BCS pairing, i.e. for IPAIRI=1 and IPAHFB=0, values of NUMBSP 
define the pairing window in which the BCS equations are solved. 



2.1. General data 
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For calculations with the HFB pairing, i.e. for IPAIRI=1 and IPAHFB=1, values of NUMBSP 
define the numbers of canonical states kept after diagonalizing density matrices. They are 
irrelevant for the pairing window, which is then defined by the cutoff energy ECUTOF. 
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2. INPUT DATA FILE 



2.2. Interaction 

This section lists keywords pertaining to definitions of coupling constants and interaction pa- 
rameters. 

2.2.1. Keyword: SKYRME-SET 



CHARACTER*4 acronym of the Skyrme force parameter set. Must start at the 13-th column of the 
data line. In version (v2.40h), valid acronyms [H] are SII, SIII, SV, SVI, SKM*, SKP, SKMP, 
SKIl, SKO, SKOP, SLY4, SLY5, SLY6, SLY7, MSKl, MSK2, MSK3, MSK4, MSK5, MSK6, 
SKX, SKXC Other sets of parameters can easily be included in the subroutine params. 



2.2.2. Keyword: SKYRME-STD 



Parameters of several standard Skyrme forces are encoded within the program. Calculation for 
a given standard force can be requested by specifying in the input data file its acronym SKYRME. 
Along with the force parameters, for each force there is encoded information on how the given 
force should be used, i.e., with which value of the parameter ^^/2m, and with which treatment 
of tensor, spin-orbit, and center-of-mass terms. For ISTAND=1, calculations are performed with 
these features set in the way specific for the given force, and the rest of the switches read on the 
same line is ignored. For ISTAND=0, the switches define nonstandard features in the following 
way: 

• For KETA_J=1 the code takes into account in the functional, and for KETA_J=0 neglects, 
the so-called tensor terms. The second option is equivalent to setting in Eq. (1-12) 
coupling constants C/=0 and Cj=Q. 

• For KETA_W=0 or 1 the code uses the traditional {Wq = Wq) or generalized (Wq ^ Wq) 
spin-orbit term, respectively, where 



For forces that use generalized spin-orbit term, option KETA_W=0 sets equal to Wq^ 
while for forces that use traditional spin-orbit, option KETA_W=1 has no effect. For 
KETA_W=2, strengths Wq and are set equal to the values of WO_INP and WOPINP, re- 
spectively, which are read under the keyword of SPIN_ORBIT. Note that if KETA_W=2 is 
used, and keyword SPIN_ORBIT is not specified in the input data file, the default values 
of WO_INP and WOPINP supersede those that are encoded within the program for the given 
Skyrme force. 



SKM* = SKYRME 



0, 1, 0, 0, = ISTAND, KETA_J, KETA_W, KETACM, KETA_M 




(4a) 
(4b) 
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For KETACM=0 the code uses the traditional one-body center-of-mass correction before 
variation, 

Ee.m.^i?f;. = -^(T), (5) 

uses the two-body center-of-mass correction after variation 



where T=—^ Hi=i is the one-body kinetic-energy operator. For KETACM=1 the code 



where -P=— «^ I^iLi is the total hnear momentum operator. Value of KETACM=2 is 
reserved for a future implementation of the two-body center-of-mass correction before 
variation. For KETACM=3 the center-of-mass correction is neglected. 

For KETA_M=0 the code uses the value of ^^/2m=20. 73620941 MeVfm^, which was encoded 
in version (vl.75r). For KETA_M=1 the code uses the value specific for the given Skyrme 
force. For KETA_M=2 the code uses the value specified in the input data file under keyword 
HBAR20VR2M. 



2.2.3. Keyword: HBAR20VR2M 

20.73620941 = HBMINP 

Value of the Ti^ /2m parameter, which is used if KETA_M=2 is set under keyword SKYRME-STD, 
and ignored otherwise. 



2.2.4. Keyword: SPIN_ORBIT 

120.0, 120.0 = WO_INP, WOPINP 

Strengths Wq and Wq of the generalized spin-orbit interaction (j4]), which are used if KETA_W=2 
is set under keyword SKYRME-STD, and ignored otherwise. 



2.2.5. Keyword: LANDAU 

000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 =LANODD, XO_LAN, X1_LAN, 

GO_LAN, GOPLAN, 
G1_LAN, GIPLAN 

Three-digit steering switch LANODD is used to specify the type of operations performed in order 
to define selected time-odd coupling constants based on values of the Landau parameters. The 
switch is followed by values of xq, xi, g^, q'q, gi, and g[, as given by equations: 

xo = C^[0]/C^[ps.t], (7a) 
xi = Cm/Cnps.t], (7b) 
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N,{2C^,[p,,,]+2Cl(3plil), (7c) 
g', = Noi2C([p.,,]+2Cf/3pli^), (7d) 
g, = -2N,C^(3plLl (7e) 



g[ = -2N,C^PpILI (7f) 




where {3 = (37rV2)2/3 and 



For LANODD=000, coupling constants are not defined by the Landau parameters, and the rest 
of parameters read on the same data hne is ignored. For LAN0DD=111, Eqs. ([7]) are solved, and 
coupling constants Co[0], Co[psat]) C*i[0], C'liPsat], Cq, and Cf are determined from the values 
of Xo, Xi, go, (7q, gi, and g[. For other values of LANODD, the following three steps are performed 
in sequence: 

1. Step one: If the rightmost digit of LANODD is equal to 1 then the coupling constants Cq 
and Cf are determined form gi and g[, Eqs. (I7e|) and (iTfl) : otherwise these two coupling 
constants are determined from the Skyrme force parameters. 

2. Step two: If the middle digit of LANODD is equal to 1 then the coupling constants CQ[psat] 
and Cf[psat] are determined form go and g'o, Eqs. (ITcl) and (l7dl) : otherwise these two 
coupling constants are determined from the Skyrme force parameters. 

3. Step three: If the leftmost digit of LANODD is equal to 1 then the coupling constants Co[0] 
and C([0] are determined form xq and xi, Eqs. fl7a|) and fiTbl) : otherwise these two coupling 
constants are determined from the Skyrme force parameters. 

All combinations of zeros and ones are allowed in LANODD. 



2.2.6. Keyword: LANDAU-SAT 

-1.0, -LO, -1.0 = HBMSAT, RHOSAT, EFFSAT 

Values of parameters /2m, psat, and m* /m, respectively, which are used when solving Eq. 
(j7]). If a negative number is read for any of these parameters, then the program uses the 
corresponding value calculated from parameters of the given Skyrme force. 



2.2.7. Keyword: EVE_SCA_TS 



1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1. 

SRHO_T, SRHO_S,SRHODT, SRHODS,SLPR_T, SLPR_S, 

STAU_T, STAU_S, 
SSCU_T, SSCU_S, 
SDIV_T, SDIV_S 



By using this item, the coupling constants corresponding to a given Skyrme parameter set can 
be arbitrarily scaled. This allows calculations with modified Skyrme functionals. The time-even 
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coupling constants in the total-sum representation, 



'-"tot 



r<P r^P 
•-"O "~ '-"1; 

2Cf, 



(9a) 
(9b) 



are multiplied by the 12 numbers from SRHO_T to SDIV_S. Variables with names ending with _T 
and _S multiply the "total" and "sum" coupling constants, respectively. Variables with names 
containing the acronyms RHO, LPR, TAU, SCU, and DIV, multiply the coupling constants with 
superscripts p, Ap, r, J, and VJ, respectively, and those with RHOD multiply the density- 
dependent part of C . Similar name convention is used for many other variables in the code 

HFODD. 



2.2.8. Keyword: ODD_SCA_TS 



Same as above but for the time-odd coupling constants. The acronyms SPI, LPS, CUR, KIS, and 
ROT correspond to coupling constants with superscripts s. As, T, j, and Vj, respectively, and 
those with SPID correspond to the density-dependent part of C"^. 



2.2.9. Keyword: EVE_SCA_PM 



Same as above but for the time-even coupling constants in the isoscalar-isovector representation. 



1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1. 

SSPI_T, SSPI_S,SSPIDT, SSPIDS,SLPS_T, SLPS_S, 



SCUR_T, SCUR_S, 
SKIS_T, SKIS_S, 
SROT_T, SROT_S 



1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1. 

SRHO_P, SRHO_M,SRHODP, SRHODM,SLPR_P, SLPR_M, 



STAU_P, STAU_M, 
SSCU_P, SSCU_M, 
SDIV_P, SDIV_M 




tot ) 



(10a) 
(10b) 



Variables with names ending with _P and _M multiply the isoscalar and isovector coupling 
constants, respectively. The total-sum scaling factors are used first, and the isoscalar-isovector 
scaling factors are used afterwards. 
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2.2.10. Keyword: ODD_SCA_PM 

1., 1., 1., 1., 1., 1, 1., 1., 1., 1., 1., 1. 

SSPI_P, SSPI_M,SSPIDP, SSPIDM,SLPS_P, SLPS_M, 

SCUR_P, SCUR_M, 
SKIS_P, SKIS_M, 
SROT_P, SROT_M 

Same as above but for the time-odd coupling constants in the isoscalar-isovector representation. 



2.2.11. Keyword: EVE_ADD_TS 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ARHO_T, ARHO_S,ARHODT, ARHODS, ALPR_T, ALPR_S, 

ATAU_T, ATAU_S, 
ASCU_T, ASCU_S, 
ADIV_T, ADIV_S 

By using this item, the couphng constants corresponding to a given Skyrme parameter set can 
be shifted by arbitrary values. The time-even coupling constants in the total-sum representation 
(1-14) are modified by adding the 12 numbers from ARHO_T to ADIV_S. By setting the scaling 
factors SRHO_T to SDIV_S equal to zero, see the keyword EVE_SCA_TS, one can input here a new 
set of the coupling constants. The name convention of variables is here the same as for many 
other variables in the code HFODD, see the keyword EVE_SCA_TS. 



2.2.12. Keyword: ODD_ADD_TS 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ASPI_T, ASPI_S,ASPIDT, ASPIDS, ALPS_T, ALPS_S, 

ACUR_T, ACUR_S, 
AKIS_T, AKIS_S, 
AROT_T, AROT_S 

Same as above but for the time-odd coupling constants. 



2.2.13. Keyword: EVE_ADD_PM 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ARHO_P, ARHO_M,ARHODP, ARHODM, ALPR_P, ALPR_M, 

ATAU_P, ATAU_M, 
ASCU_P, ASCU_M, 
ADIV_P, ADIV_M 

Same as above but for the time-even coupling constants in the isoscalar-isovector representation 
( JTO|) . The total-sum additive factors are used first, and the isoscalar-isovector additive factors 
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are used afterwards. 



2.2.14. Keyword: ODD_ADD_PM 

0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. 

ASPI_P, ASPI_M,ASPIDP, ASPIDM, ALPS_P, ALPS_M, 

ACUR_P, ACUR_M, 
AKIS_P, AKIS_M, 
AROT_P, AROT_M 

Same as above but for the time-odd coupling constants in tlie isoscalar-isovector representation 



2.2.15. Keyword: G_SCALING 

1.0, 1.0 = FACTGN, FACTGP 

For IPABCS=1 the code HFODD solves the BCS equations with the neutron and proton pairing 
strengths defined in Ref. [5]. These values can be modified by defining here the multiplicative 
factors FACTGN and FACTGP for neutrons and protons, respectively. 



2.2.16. Keyword: PAIR_MATRI 

1, 0, 0, = IDESTA, IDEMID, IDESTO, IDEDIS 
For IDESTA=1, the pairing matrix elements, required for the BCS pairing calculations with 
state-dependent pairing gaps (IPABCS=3), are calculated in the first iteration. In the present 
version (v2.40h), only the value of IDESTA=1 is allowed, because the option of storing the 
pairing matrix elements is not yet implemented. For IDEMID=1, the pairing matrix elements 
are calculated in the middle iterations (all between the first and the last one), for IDEST0=1 in 
the last iteration, and/or for IDEDIS=n, in every n-th iteration, counting from the first one. 



2.2.17. Keyword: INI_FERMI 

-8.0, -8.0 = FERINI(O), FERINI(l) 



2.2.18. Keyword: INI_DELTA 

1.0, 1.0 = DELINI(O), DELINI(l) 

For ICONTI=0 or IC0NTI=1 and IPCONT=0, the HFB calculations are started by using the initial 
values of the neutron and proton Fermi energies, FERINI(O) and FERINI(l), and pairing gaps, 
DELINI(O) and DELINI(l). For IC0NTI=1, the initial values of pairing gaps overwrite the 
HFB pairing potentials with constant values of DELINI(O) and DELINI(l) for neutrons and 
protons, respectively. These constant potentials are ignored after the first iteration, i.e., in the 
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first iteration, the mixing of previous and current potentials (see keyword SLOW-PAIR) is not 
performed. A possibility of restarting calculations with nonzero pairing is very useful in case 
the pairing would have vanished in a former run. For ICONTI=0 or IPCONT=0, values of FERINI 
and DELINI are ignored. 



2.2.19. Keyword: FIXDELTA_N 

1.0, = DELFIN, IDEFIN 

For IDEFIN=1, the HFB pairing calculations (for IPAHFB=1) or BCS pairing calculations (for 
IPABCS=2) are performed with a fixed value of the neutron pairing gap equal to DELFIN. For 
IDEFIN=0, value of DELFIN is ignored and the zero-range density-dependent pairing force f ill I) 
is used for neutrons. IDEFIN=1 requires IPAHFB=1 or IPABCS=2. 



2.2.20. Keyword: FIXDELTA_P 

1.0, = DELFIP, IDEFIP 
Same as above but for the proton pairing gap. 



2.2.21. Keyword: PAIRNFORCE 

-200.0, 0.0, 1.0 = PRHO_N, PRHODN, POWERN 

Parameters Vq, Vi, and a, respectively, of the zero-range density- dependent pairing force defined 
by the pairing form-factor: 



which is used in the HFB pairing or BCS pairing calculations for neutrons. In case the values 
of Vi and a allow it, the code calculates the value of the saturation density po that gives the 
equivalent form of the pairing form-factor ( ITTi) . In case po cannot be calculated, the codes set 
its value to 1; po is calculated only for the purpose of information, while internally the code 
uses only the value of Vi. For IDEFIN=1, or for IPAHFBt^I and IPABCSt^S, parameters PRHO_N, 
PRHODN, and POWERN are ignored. 




(11) 



2.2.22. Keyword: PAIRPFORCE 

-200.0, 0.0, 1.0 
Same as above but for the proton pairing force. 



= PRHO_P, PRHODP, POWERP 
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2.2.23. Keyword: PAIR_FORCE 

-200.0, 0.0, 1.0 = PRHO_T, PRHODT, POWERT 

Same as above but for the neutron and proton pairing forces. Tliis keyword is equivalent to 
using the two keywords simultaneously, PAIRNFORCE and PAIRPFORCE, with identical parameters 
for neutrons and protons. 



2.2.24. Keyword: PAIRNINTER 

-200.0, 0.16, 1.0 = PRHO_N, PRHOSN, POWERN 

Parameters Vq, po, and a, respectively, of the zero-range density-dependent pairing force (fTTl) . 
which is used in the HFB pairing or BCS pairing calculations for neutrons. In case the values 
of Po and a allow it, the code calculates the value of Vi that gives the equivalent form of the 
pairing form-factor ffTTl) . In case Vi cannot be calculated, the code stops. For IDEFIN=1, or for 
IPAHFBt^I and IPABCSt^S, parameters PRHO_N, PRHOSN, and POWERN are ignored. 



2.2.25. Keyword: PAIRPINTER 

-200.0, 0.16, 1.0 = PRHO_P, PRHOSP, POWERP 
Same as above but for the proton pairing force. 



2.2.26. Keyword: PAIR_INTER 

-200.0, 0.16, 1.0 = PRHO_T, PRHOST, POWERT 

Same as above but for the neutron and proton pairing forces. This keyword is equivalent to 
using the two keywords simultaneously, PAIRNINTER and PAIRPINTER, with identical parameters 
for neutrons and protons. 



2.2.27. Keyword: CUTOFF 

60.0 = ECUTOF 

The cutoff energy Cmax for summing up contributions of the HFB quasiparticle states to density 
matrices, see Sec. IV-2.5. Ignored for IPAHFB=0. 



2.2.28. Keyword: YUKAWATERM 

0.7045, 4.7565, 1.0, 0.0, 0.0, 0.0, 1, 

PIMASS, PNMASS, YUKAGT, YUKAGO, YUKAGl, YUKAG2, 

lYUTYP, I_YUKA 

For I_YUKA>0, the code calculates the average values of the time-reversal- and parity-violating 
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Yukawa interaction (VI-59), with the pion mass (m^) of PIMASS and the nucleon mass (mjv) of 
PNMASS. If values of zero are read, variables PIMASS and PNMASS remain unchanged. Variables 
YUKAGT, YUKAGO, YUKAGl, and YUKAG2 correspond, respectively, to the couphng constants 
go, cji, and g2, see Eq. (VI-59). For I_YUKA=2 or 3, the direct matrix elements of the Yukawa 
interaction (VI-59) are, in addition, added to the self-consistent mean field. For I_YUKA=2 or 
4, the exchange matrix elements are added. For IYUTYP=1, expression (VI-61) is used, while 
for IYUTYP=2, an analogous six-Gaussian expression is used without the short-range correction 
(VI-58), that is, for /(r) = 1. For I_YUKA=0, all these input data are ignored and the Yukawa 
interaction is not taken into account. 

Figures [U and [2] show the actual CPU times required in calculations using the Skyrme 
and Yukawa interactions, respectively. At A'o = 10, the latter take more than two orders of 
magnitude longer and, moreover, their CPU times scale like Nq rather than A^q . A local EDF 
clearly makes for easier computing. 




Figure 1: The HFODD CPU times required for calculations that use the standard Skyrme EDF, 
shown as a function of the number of HO shells Nq. The doubly logarithmic scale in the Figure, 
shows that these times scale as An. 



2.2.29. Keyword: COULOMBPAR 

7, 1, 1 = ICOTYP, ICOUDI, ICOUEX 
For ICOUDI=0, 1, or 2, the Coulomb direct energy and Coulomb direct mean field are neglected, 
calculated by using the Green- function method, see Section 1-5, or calculated by using the 
Gaussian-expansion method, see Section VI-2.10, respectively. Similarly, for ICOUEX=0, 1, or 
2, the Coulomb exchange energy and Coulomb exchange mean field are neglected, calculated 
by using the Slater approximation (1-19) or calculated exactly by using the Gaussian-expansion 
method, see Section VI-2.10, respectively. For the Gaussian-expansion method, that is, for 
IC0UDI=2 or IC0UEX=2, positive values of ICOTYP denote the numbers of Gauss-Legendre 
nodes used in the integral of Eq. (VI-103). For IC0UDI=2 or IC0UEX=2, the iteration can later 
be smoothly continued (IFC0NT=1, see the keyword CONTFIELDS) only by saving the matrix 
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Figure 2: Same as in Figured] but for calculations that use the Yukawa interaction. The CPU 



times scale as Nq. 



elements of the mean field, that is, by requesting IWRIFI=1, see the keyword FIELD_SAVE. 
Therefore, IC0UDI=2 or IC0UEX=2 and IC0NTI=1 requires IFC0NT=1. 

In Fig. [Hlthe error in the exact Coulomb exchange energy is plotted as function of the number 
of Gauss-Legendre nodes Nq. A quite precise estimate is obtained for Nc = 7 Gauss- Legendre 
nodes and the machine precision is obtained by doubling this number (along with doubling the 
CPU time). For Nc = 7, these CPU times are shown in Fig. |H which illustrates the fact that 
these times scale with the number of HO shells as N^. 
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Figure 3: Errors in the exact Coulomb exchange energy plotted in function of the number of 
Gauss-Legendre nodes Nc- 
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No 

Figure 4: Same as in Figure [1] but for the calculations that determine the exact Coulomb 
exchange energies. The CPU times scale as N^. 



2.2.30. Keyword: INSERT_HO 

= IPOTHO 

For IP0TH0=1, an external HO potential is added to the self-consistent mean field. Parameters 
of the potential are identical to those defining the HO basis. 



2.3. Symmetries 



25 



2.3. Symmetries 

This section lists keywords pertaining to selection of conserved, broken, and restored symme- 
tries. 



2.3.1. Keyword: SIMPLEXY 

1 = ISIMPY 

For ISIMPY=1, calculation with y-simplex conserved are performed, while for ISIMPY=0 the 
simplex is broken. Value of ISIMPY must be consistent with switches IPARTY, ISIGNY, ISIMTY, 
and IROTAT, see Tables □ and El 



2.3.2. Keyword: SIGNATUREY 

1 = ISIGNY 

For ISIGNY=1, calculation with signature conserved are performed, while for ISIGNY=0 the 
signature is broken. Value of ISIGNY must be consistent with switches ISIMPY, IPARTY, ISIMTX, 
and ISIMTZ, see Tables □ and O 



2.3.3. Keyword: PARITY 

-1 = IPARTY 

For IPARTY=1, calculation with parity conserved are performed, while for IPARTY=0 the parity 
is broken. Value of IPARTY must be consistent with switches ISIMPY and ISIGNY, see Table 
[H For IPARTY=— 1 (the compatibility mode), the code sets IPARTY=ISIGNY and requires that 
ISIMPY=1 



2.3.4. Keyword: ROTATION 

1 = IROTAT 

For IR0TAT=1, calculation with time-reversal breaking are performed, while for IROTAT=0 the 
time-reversal symmetry is conserved. In the latter case the calculations are performed only 
for one value of the simplex, s=+i, which gives almost twice shorter CPU times. IROTAT=0 
is incompatible with providing a non-zero value of the angular frequency or with attempting 
calculation for odd or odd-odd nuclei without pairing, IPAIRI=0. IR0TAT=1 is incompatible 
with IPABCS>0, i.e., the rotating solutions can be obtained for the HFB pairing but not for the 
BCS pairing. 
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Table 1: Primary set of conserved and nonconserved symmetries allowed in the code hfodd 
version (v2.40h). 



Option 


Symmetries 


Switches 


Sy 


Ry 


P 


ISIMPY 


ISIGNY 


IPARTY 


PI 


conserved 


conserved 


conserved 


1 


1 


1 


P2 


conserved 


nonconserved 


nonconserved 


1 








P3 


nonconserved 


conserved 


nonconserved 





1 





P4 


nonconserved 


nonconserved 


conserved 








1 


P5 


nonconserved 


nonconserved 


nonconserved 












Table 2: Secondary set of conserved and nonconserved symmetries allowed in the code hfodd 
version (v2.40h). 



Option 


Symmetries 


Switches 


Sy 


f 


^y 


ISIMPY 


ITIREV 


ISIMTY 


SI 


conserved 


conserved 


conserved 


1 


1 


1 


S2 


conserved 


nonconserved 


nonconserved 


1 








S3 


nonconserved 


conserved 


nonconserved 





1 





S4 


nonconserved 


nonconserved 


conserved 








1 


S5 


nonconserved 


nonconserved 


nonconserved 












Table 3: Ternary set of conserved and nonconserved symmetries allowed in the code hfodd 
version (v2.40h). 



Option 


Symmetries 


Switches 


Ry 


si 


qT 


ISIGNY 


ISIMTX 


ISIMTZ 


Tl 


conserved 


conserved 


conserved 


1 


1 


1 


T2 


conserved 


nonconserved 


nonconserved 


1 








T3 


nonconserved 


conserved 


nonconserved 





1 





T4 


nonconserved 


nonconserved 


conserved 








1 


T5 


nonconserved 


nonconserved 


nonconserved 
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2.3.5. Keyword: TIMEREVERS 

= ITIREV 

For ITIREV=1, calculation with time-reversal conserved are performed, while for ITIREV=0 this 
symmetry is broken. This switch is used only as a convenient replacement of switch IROTAT; 
ITIREV=1 is equivalent to IROTAT=0 and ITIREV=0 is equivalent to IR0TAT=1. 



2.3.6. Keyword: TSIMPLEX_Y 

-1 = ISIMTY 

For ISIMTY=1, calculation with y-simplex-^ conserved are performed, while for ISIMTY=0 this 
symmetry is broken. Value of ISIMTY must be consistent with switches ISIMPY and IROTAT, 
see Table H For ISIMTY=-1 (the compatibility mode), the code sets ISIMTY to 1-IROTAT and 
requires that ISIMPY=1. 



2.3.7. Keyword: TSIMPLEXES 

1, 1 = ISIMTX, ISIMTZ 

For ISIMTX=1 and/or ISIMTZ=1, calculation with conserved symmetries given by x-simplex^ 
and/or z-simplex"^ are performed, respectively. These symmetries are broken for ISIMTX=0 and 
ISIMTZ=0. Values of ISIMTX and ISIMTZ must be consistent with that of ISIGNY, see Table H 



2.3.8. Keyword: TSIMPLEX3D 

1, -1, 1 = ISIMTX, ISIMTY, ISIMTZ 

This keyword allows to simultaneously input all the three switches that define the three T- 
simplexes. It is equivalent to using keywords TSIMPLEXES and TSIMPLEX_Y together. 



2.3.9. Keyword: PAIRING 

= IPAIRI 

For IPAIRI=0, calculation without pairing correlations (pure HF) are performed, while for 
IPAIRI=1 the pairing correlations are included. In the latter case, for IPAHFB=0 or 1 pairing is 
included within the BCS or HFB method, respectively. IPAIRI=1 requires either IPAHFB=1 or 
IPABCS>0 and is incompatible with simultaneously setting IPAHFB=1 and IPABCS>0, i.e., the 
HFB pairing and BCS pairing cannot be used simultaneously. IPAIRI=1 is incompatible with 
IR0TAT=1, unless IPAHFB=1, i.e., the rotating solutions can be obtained for the HFB pairing 
but not for the BCS pairing. 
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2.3.10. Keyword: HFB 

= IPAHFB 

For IPAHFB=1, calculation with pairing correlations included by solving the HFB equation are 
performed. IPAHFB=1 requires IPAIRI=1. 



2.3.11. Keyword: BCS 

-1 = IPABCS 

Parameter IPABCS defines the type of BCS pairing calculations. For IPABCS=0, 1, 2, or 3, the 
BCS pairing calculations are, respectively, not performed, performed with the seniority pairing 
force, performed with fixed pairing gaps, or performed with the state-dependent pairing gaps. 
Value of IPABCS=— 1 is allowed for the sake of compatibility with earlier versions of the code, 
before (v2.13f). Then, IPABCS=0 is set for IPAHFB=1 and IPABCS=IPAIRI is set for IPAHFB=0. 

For IPABCS=1 the code HFODD solves the BCS equations with the neutron and proton pairing 
strengths defined in Ref . [9] . These values can be modified by defining the multiplicative factors 
FACTGN and FACTGP for neutrons and protons, respectively. For IPABCS=2, the BCS pairing 
calculations are performed with fixed values of the neutron or proton pairing gaps equal to 
DELFIN or DELFIP, respectively. For IPABCS=3, the BCS pairing calculations are performed 
with state-dependent pairing gaps corresponding to the pairing matrix elements (see the key- 
word PAIR_MATRI) calculated for the contact forces defined in the same way as for the HFB 
pairing calculations. 

Positive values of IPABCS require IPAIRI=1 and IPAHFB=0 and are incompatible with rotations 
IR0TAT=1 or broken simplex ISIMPY=0. IPABCS=2 requires IDEFIN=IDEFIP=1. 



2.3.12. Keyword: FILSIG_NEU 

2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, = 

KPFILG (0,0,0), KPFILG (0,1,0), KPFILG(1 , , 0) , KPFILG(1 , 1 , 0) , 
KHFILG (0,0,0), KHFILG (0,1,0), KHFILG(1 , , 0) , KHFILG(1 , 1 , 0) , 
KOFILG(0,0,0), KOFILG (0,1,0), K0FILG(1 , , 0) , K0FILG(1 , 1 ,0) 

These parameters govern the filling approximation for the neutron single-particle parity-signa- 
ture configurations. Matrices KPFILG contain the indices of particle states in the four parity- 
signature blocks denoted by (+,+), (+,—), (— 5+), and (—,—), of given (paritysignature) com- 
binations, i.e., (7r,r) = and respectively. Matrices KHFILG 
contain analogous indices of hole states, and matrices KOFILG contain numbers of particles put 
into the states between KHFILG and KPFILG, by using for them partial occupation factors of 
K0FILG/(KPFILG-KHFILG-hl)/2. For KOFILG=0, in the given parity-signature block the filhng 
approximation is inactive. The filling approximation is incompatible with pairing correlations, 
IPAIRI=1. 



2.3. Symmetries 
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2.3.13. Keyword: FILSIG_PRO 

2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, = 

KPFILG(0,0,1), KPFILG(0,1,1), KPFILG(1 , , 1) , KPFILG(1 , 1 , 1) , 
KHFILG(0,0,1), KHFILG(0,1,1), KHFILG(1 , , 1) , KHFILG(1 , 1 , 1) , 
KQFILG(0,0,1), K0FILG(0,1,1), K0FILG(1 , , 1) , K0FILG(1 , 1 , 1) 

Same as above but for the proton single-particle parity-signature configurations. 



2.3.14. Keyword: INI_INVERS 

0, = INIINV, INIKAR 

Allowed values of INIKAR=0, 1, 2, or 3 and INIINV=0, 1, 2, or 3 correspond to the following 
initial Djy, transformations: 





INIKAR=0 INIKAR=1 INIKAR=2 INIKAR=3 


INIINV=0 
INIINV=1 
INIINV=2 
INIINV=3 


/ Rx Ry Rz 

P Sx Sy Sz 

T -^x Rl 
pT cT qT qT 

r 



where P is the space inversion, T is the time-reversal symmetry, Rk is the signature (rotation 
by angle vr about the k = x,y, or z axis), and P^ = PT, Sk = PRk (simplex), i?^ = TRk 
(/c-signature-^), and S^ = TSk {k-simplex^). For INIKAR=INIINV=0, no transformation is per- 
formed and this option is inactive. 

Transformations are performed at the level of the densities, after the first iteration. As a secu- 
rity measure, nonzero values of INIKAR and INIINV require N0ITER=1. Such values also require 
SLOWEV=SLOWOD=SLOWPA=0.0, so as not to mix the old and new potentials corresponding to 
the original and transformed densities, respectively. They are also incompatible with IPRGCMt^O. 

A given initial Djj^ transformation must be accompanied by the correspondingly broken sym- 
metries, that is, INIINV=1 or 3 requires IR0TAT=1 and INIINV=2 or 3 requires IPARTY=0. 



2.3.15. Keyword: INI_ROTAT 

0.0, 0.0, 0.0, = ALPINI, BETINI, GAMINI, INIROT 
For INIR0T=1, the wave functions are rotated by the Euler angles a, (3, and 7 corresponding, 
respectively, to ALPINI, BETINI, and GAMINI (all in degrees). Transformations are performed 
at the level of the densities, after the first iteration. As a security measure, INIR0T=1 requires 
N0ITER=1. It also requires SLOWEV=SLOWOD=SLOWPA=0.0, so as not to mix the old and new 
potentials corresponding to the original and transformed densities, respectively. INIR0T=1 is 
incompatible with IPRGCMt^O and requires a spherical HO basis of hujx = huy = Tiujz. 
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A rotation about the z axis must be accompanied by the broken signature, that is, ALPINIt^O 
or GAMINIt^O requires ISIGNY=0. 



2.3.16. Keyword: PROJECTGCM 

0,0,0, 1,1,0, 1,1,0 

IPRGCM, IPROMI, IPROMA, NUAKNO, NUBKNO, KPROJE, 

IFRWAV, ITOWAV, IWRWAV 

For IPRGCM=1 and 2, the code calculates the diagonal and non-diagonal GCM kernels, respec- 
tively. In addition, for NUAKNOt^I or NUBKNOt^I, the AMP kernels are calculated, as described 
in Section VI-2.1, although in the present version (v2.40h) this option is not yet available for 
non-diagonal kernels of IPRGCM=2. The angular momentum projection (AMP) is performed 
for doubled angular momenta from IPROMI to IPROMA and requires a spherical HO basis of 
Uu)^ = hujy = huz- For even (odd) particle number IN_FIX+IZ_FIX, IPROMI, IPROMA, and 
KPROJE must be even (odd). For IPRGCM=1 and 2, the present version (v2.40h) also requires 
IR0TAT=1 and IPAIRI=0. 

The AMP has been tested up to the values of angular momenta of 70^, and therefore, IPROMA 
must not be larger than 2*70=140. After further tests, higher values could be used by increas- 
ing the parameter JMAX=70 in function DSMALG, which calculates the Wigner functions, see 
Eqs. (VI-1), (VI-4), and (VI-5). 

NUAKNO is the number of Gauss-Tchebyschev nodes, which are used to perform integrations over 
the a and 7 Euler angles. For NUAKN0=1, these integrations are not performed (ID AMP) and 
the states are assumed to be axial with the doubled projection of the angular momentum on 
the z axis equal to KPROJE. NUBKNO is the number of Gauss-Legendre nodes, which are used to 
perform integrations over the (3 Euler angle. For NUAKN0>1 and NUBKN0>1, a full 3D AMP is 
performed and the value of KPROJE is ignored. NUAKN0>1 requires ISIMPY=0 and ISIGNY=0. 
IPROMA must be larger than the absolute value of KPROJE. 

For IPRGCM=2, the code calculates the GCM kernels between the states labeled by three-digit 
indices from "000" to "999". Indices of the "left" states vary between IFRWAV and ITOWAV, 
and these states are read from the disc. The index of the "right" state equals to ITOWAV, and 
this state is equal to the current state. In addition, for IWRWAV=1, the current state is saved 
on the disc with the index of ITOWAV. This allows for a simultaneous calculation of the "right" 
state along with calculating its kernels with all previously calculated and stored "left" states. 
The states can also be stored on disc without calculating kernels in the given run, that is, by 
setting the value of ITOWAV along with IWRWAV=1 and IPRGCM=0. Names of files on the disc 
are composed by concatenating the three-digit index, "-", and FILWAV. 



2.3.17. Keyword: SAVEKERNEL 

= ISAKER 



2.3. Symmetries 
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For ISAKER=1, the code attempts reading tlie kernel files Nxxx-Lyyy-Rzzz-/ /FILKER, where / / 
denotes concatenated strings. The three-digit indices are: 

• XXX is the consecutive index of the kernel file, 

• yyy is the index of the left wave function, 

• zzz is the index of the right wave function. 

In the work directory, the file names for all indices xxx are scanned, starting from 001. The 
kernels stored in these files are read into memory and are not recalculated. The kernels that 
have not been found in the kernel files are calculated and stored in the kernel file with the 
lowest available index xxx. In this way, one can submit many parallel jobs, see the keyword 
PARAKERNEL, that calculate kernels for different values of the Euler angles a, and 7. The 
results are then collected in different kernel files with indices xxx attributed automatically. If 
any of the jobs is terminated before completing its task, the same input data can be resubmitted 
and the calculation automatically continues from the point where it has been interrupted. Once 
all the kernels will have been calculated, which requires a large CPU time, the AMP can be 
performed within a very small CPU time by reading, again automatically, all the created kernel 
files. ISAKER=1 requires IPRGCM>0 



2.3.18. Keyword: CHECKERNEL 

1 = ICHKER 

The names of kernel files are saved within these files. As a security measure, when reading the 
kernel files, their names are cross-checked against the saved information. This cross-checking 
can be deactivated by using ICHKER=0. This option is useful whenever the kernel files have 
been renamed for any reason. 



2.3.19. Keyword: PARAKERNEL 

0, 1, 1, 1, 1 = IPAKER, NUASTA, NUASTO, NUGSTA, NUGSTO 
For IPAKER=1, the code only calculates kernels for different values of the Euler angles a, (3, and 
7 and the AMP is suspended. Calculations are performed for the Gauss-Tchebyschev nodes 
in the Euler angle a from NUASTA to NUASTO, for those in the Euler angle 7 from NUGSTA to 
NUGSTO, and for all the Gauss- Legendre nodes in the Euler angle f3, that is, from 1 to NUBKNO. 
For IPAKER=1, to save memory, the code can be compiled with IPARAL=1. IPAKER=1 requires 
IPRGCM>0 and ISAKER=1. Values of NUASTA, NUASTO, NUGSTA, NUGSTO must all be between 1 
and NUAKNO and must be ordered as NUASTA<NUASTO and NUGSTA<NUGSTO. 



2.3.20. Keyword: TRANSITION 

2, 1, = NMURED, NMARED, NSIRED 
Maximum numbers of transition electric, transition magnetic, and transition surface or Schiff 
moments, respectively, for which proton kernels and reduced matrix elements are calculated, 
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printed, and stored in the kernel files. NMARED and NSIRED must not exceed NMURED. For 
NMURED=0, NMARED=0, or NSIRED, the corresponding proton kernels and reduced matrix ele- 
ments are not calculated. 



2.3.21. Keyword: CUTOVERLAP 

0, 10"^°, 1. = ICUTOV, CUTOVE, CUTOVF 
For ICUTOV=0, parameters CUTOVE and CUTOVF are ignored and the collective states for the 
iC-mixing calculation, Eq. (VI-10), are selected by their norm eigenvalues being larger than 
the negative of the smallest norm eigenvalue OVEMIN. For ICUT0V=1, the collective states are 
selected by their norm eigenvalues being larger than CUTOVE+CUTOVF*OVEMIN. 



2.3.22. Keyword: LIPKIN 

0, = LIPKIN, LIPKIP 

For LIPKIN=1 and/or LIPKIP=1, the Lipkin-Nogami corrections are included for neutrons 
and/or protons, respectively, see Section VI-29. In the present version (v2.40h), LIPKIN=1 or 
LIPKIP=1 requires IPAHFB=1. 



2.3.23. Keyword: INI_LIPKIN 

0.1, 0.1 = FE2INI(0), FE2INI(1) 
For IPCONT=0, or IPC0NT=1 and ILCONT=0, the Lipkin-Nogami calculations are started by 
using the initial values of the neutron and proton Lipkin-Nogami parameters A2 (VI-98), 
FE2INI(0) and FE2INI(1), respectively. For IPC0NT=1 and ILC0NT=1, values of FE2INI are 
ignored. 



2.3.24. Keyword: FIXLAMB2_N 

0.1, = FE2FIN, IF2FIN 
For IF2FIN=1, the neutron Lipkin-Nogami calculations are performed by using the fixed value 
of the neutron Lipkin-Nogami parameter A2 (VI-98) equal to FE2FIN. IF2FIN=1 requires 
LIPKIN=1. 



2.3.25. Keyword: FIXLAMB2_P 

0.1, = FE2FIP, IF2FIP 
Same as above but for the proton Lipkin-Nogami calculations. 



2.3. Symmetries 
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2.3.26. Keyword: SLOWLIPKIN 

0.5 = SLOWLI 

SLOWLI gives tlie value of tlie mixing fraction e used for the Lipkin-Nogami parameter A2 (VI- 
98), in analogy to the SLOWEV, SLOWOD and SLOWPA parameters. 



2.3.27. Keyword: FIXFERMI_N 

-8.0, = FERFIN, IFEFIN 
For IFEFIN=1, the HFB pairing calculations are performed with a fixed value of the neutron 
Fermi energy equal to FERFIN. In the present version (v2.40h), IFEFIN=1 requires IPAHFB=1. 



2.3.28. Keyword: FIXFERMI_P 

-8.0, = FERFIP, IFEFIP 

Same as above but for the proton HFB pairing calculations. 
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2.4. Configurations with no conserved symmetries 

This section lists keywords pertaining to definitions of configurations in the case when no 
symmetries are conserved, and the single-particle or single-quasiparticle space is divided into 
two charge blocks (neutrons and protons). 



2.4.1. Keyword: PHNONE_NEU 

1, 00, 00 = NUPAHO, KPNONE, KHNONE 

Neutron particle-hole excitations pertaining to the situation when all neutrons are in one com- 
mon block (no simplex, signature, or parity is conserved). NUPAHO is the consecutive number 
from 1 to 5 (up to five sets of excitations can be specified in separate items). Particles are re- 
moved from the KHNONE-th state and put in the KPNONE-th state. At every stage of constructing 
excitations the Pauli exclusion principle has to be respected (particle removed from an occupied 
state and put in an empty state). Values equal zero have no effect. Note that for all neutrons 
sitting in one common block the reference configuration from which the particle-hole excita- 
tions are counted is defined by the total number of neutrons. These particle-hole excitations 
are ignored unless ISIMPY=0, IPARTY=0, and IPAIRI=0. 

2.4.2. Keyword: PHNONE_PRO 

1, 00, 00 = NUPAHO, KPNONE, KHNONE 

Same as above but for the proton particle-hole excitations. 

2.4.3. Keyword: DIANON_NEU 

2, 1, = KPFLIZ, KHFLIZ, KOFLIZ 

The diabatic blocking of neutron single-particle configurations in the situation when all neutrons 
are in one common block. The rules to define the blocking are analogous to those described for 
the conserved parity, keyword DIAPAR_NEU. 



2.4.4. Keyword: DIANON_PRO 

2, 1, = KPFLIZ, KHFLIZ, KOFLIZ 

Same as above but for the diabatic blocking of proton single-particle configurations in the 
situation when all protons are in one common block. 



2.4.5. Keyword: BLOCKSIZ_N 

1, = INSIZN, IDSIZN 



2.4. Configurations with no conserved symmetries 
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For |IDSIZN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when no symmetries are conserved, see Section VI-2.7. For IDSIZN=+1 or —1, the blocked 
quasiparticle state is selected by having the largest overlap with the INSIZN-th neutron single- 
particle eigenstate of the HFB mean-field Routhian or with its time-reversed partner, respec- 
tively. Note that for rotating states, the time-reversed eigenstate is not necessarily an eigenstate 
of the Routhian. |IDSIZN|=1 requires ISIMPY=0, IPARTY=0, IPAHFB=1, and IR0TAT=1. 



2.4.6. Keyword: BLOCKSIZ_P 

1, = INSIZP, IDSIZP 

Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 
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2.5. Configurations with conserved simplex 

This section lists keywords pertaining to definitions of configurations in the case when simplex 
is conserved, and the single-particle or single-quasiparticle space is divided into four charge- 
simplex blocks (neutrons and protons with simplexes s=±i). 



2.5.1. Keyword: VACSIM_NEU 

43, 43 = KVASIM(O), KVASIM(l) 

Numbers of lowest neutron states occupied in the two simplex blocks, denoted by (+) and (— ), 
of given simplexes, s=+i and s=—i, respectively. These numbers define the simplex reference 
configuration from which the particle-hole excitations are counted. The definitions of simplex 
reference configuration and excitations are ignored unless ISIMPY=1, ISIGNY=0, and IPAIRI=0. 



2.5.2. Keyword: VACSIM_PRO 

33, 33 = KVASIM(O), KVASIM(l) 

Same as above but for the numbers of proton states. 

2.5.3. Keyword: PHSIMP_NEU 

1, 00, 00, 00, 00 = NUPAHO, KPSIMP, KPSIMM, KHSIMP, KHSIMM 

Neutron particle-hole excitations in the simplex blocks. NUPAHO is the consecutive number from 
1 to 5 (up to five sets of excitations can be specified in separate items). Particles are removed 
from the KHSIMP-th state in the (+) block and from the KHSIMM-th state in the (— ) block, 
and put in the KPSIMP-th state in the (-I-) block and in the KPSIMM-th state in the (— ) block. 
At every stage of constructing excitations the Pauli exclusion principle has to be respected 
(particle removed from an occupied state and put in an empty state). Values equal zero have 
no effect. In practice, reasonable excitations can only be constructed by consulting the printed 
lists of single-particle states with their consecutive numbers in given blocks. These particle-hole 
excitations are ignored unless ISIMPY=1, ISIGNY=0, and IPAIRI=0. 



2.5.4. Keyword: PHSIMP_PRO 

1, 00, 00, 00, 00 = NUPAHO, KPSIMP, KPSIMM, KHSIMP, KHSIMM 
Same as above but for the proton particle-hole excitations. 



2.5. Configurations with conserved simplex 
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2.5.5. Keyword: DIASIM_NEU 

2,2, 1,1, 0,0= KPFLIM(0,0), KPFLIM(1,0), 

KHFLIM(0,0), KHFLIM(1,0), 
KOFLIM(0,0), K0FLIM(1,0) 

The diabatic blocking of neutron single-particle simplex configurations. Matrices KPFLIM con- 
tain the indices of the particle states in the two simplex blocks denoted by (+) and (— ), of given 
simplex values, i.e., s=+i and —i, respectively. Matrices KHFLIM contain analogous indices of 
the hole states, and matrices KOFLIM define types of blocking in analogy to those shown in 
Table H 



2.5.6. Keyword: DIASIM_PRO 

2,2, 1,1, 0,0= KPFLIM(0,1), KPFLIM (1,1), 

KHFLIM(0,1), KHFLIM(1,1), 
K0FLIM(0,1), K0FLIM(1,1) 

Same as above but for the diabatic blocking of proton single-particle simplex configurations. 



2.5.7. Keyword: BLOCKSIM_N 

1, 0, = INSIMN, IRSIMN, IDSIMN 
For |IDSIMN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when simplex is conserved, see Section VI-2.7. For IDSIMN=-|-1 or —1, the blocked quasiparticle 
state is selected by having the largest overlap with the INSIMN-th neutron single-particle eigen- 
state of the HFB mean-field Routhian in a given simplex block or with its time-reversed partner, 
respectively. The simplex of the state, +i or —i, is defined by IRSIMN=0 or 1, respectively. 
I IDSIMN] =1 requires ISIMPY=1, IPARTY=0, IPAHFB=1, and IR0TAT=1. 



2.5.8. Keyword: BLOCKSIM_P 

1, 0, = INSIMP, IRSIMP, IDSIMP 
Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 



38 



2. INPUT DATA FILE 



2.6. Configurations with conserved parity 

This section lists keywords pertaining to definitions of configurations in the case when parity is 
conserved, and the single-particle or single-quasiparticle space is divided into four charge-parity 
blocks (neutrons and protons with parities 7r=±l). 



2.6.1. Keyword: VACPAR_NEU 

44, 42 = KVASIQ(O), KVASIQ(l) 

Numbers of lowest neutron states occupied in the two parity blocks, denoted by (+1) and 
(— 1), of given parities, 7r=+l and 7r=— 1, respectively. These numbers define the parity ref- 
erence configuration from which the particle-hole excitations are counted. The definitions of 
parity reference configuration and excitations are ignored unless IPARTY=1, or IPARTY=— 1 and 
ISIMPY=ISIGNY=1. They are also ignored for ISIMPY=1 and IPAIRI=1. 



2.6.2. Keyword: VACPAR_PRO 

32, 34 = KVASIQ(O), KVASiq(l) 
Same as above but for the numbers of proton states. 



2.6.3. Keyword: PHPARI_NEU 

1, 00, 00, 00, 00 = NUPAHO, KPSIQP, KPSIQM, KHSIQP, KHSIQM 

Neutron particle-hole excitations in the parity blocks. NUPAHO is the consecutive number from 
1 to 5 (up to five sets of excitations can be specified in separate items). Particles are removed 
from the KHSIQP-th state in the (+1) block and from the KHSIQM-th state in the (—1) block, and 
put in the KPSIQP-th state in the (+1) block and in the KPSIQM-th state in the (-1) block. At 
every stage of constructing excitations the Pauli exclusion principle has to be respected (particle 
removed from an occupied state and put in an empty state). Values equal zero have no effect. 
These particle-hole excitations are ignored unless ISIMPY=0, IPARTY=1, and IPAIRI=0. 



2.6.4. Keyword: PHPARI_PRO 

1, 00, 00, 00, 00 = NUPAHO, KPSIQP, KPSIQM, KHSIQP, KHSIQM 

Same as above but for the proton particle-hole excitations. 



2.6. Configurations with conserved parity 
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2.6.5. Keyword: DIAPAR_NEU 

2,2, 1,1, 0,0= KPFLIQ(0,0), KPFLIQ(1,0), 

KHFLIQ(0,0), KHFLIQ(1,0), 
KOFLIQ(0,0), K0FLIQ(1,0) 

The diabatic blocking of neutron single-particle parity configurations. Matrices KPFLIQ contain 
the indices of the particle states in the two parity blocks denoted by (+1) and (—1), of given 
parities, i.e., 7r=+l and —1, respectively. Matrices KHFLIQ contain analogous indices of the 
hole states, and matrices KOFLIQ define the type of blocking according to definitions shown in 
Table H 



Table 4: Values of parameter KOFLIQ that specify type of blocking in case of conserved parity. 
Analogous values define type of blocking for other conserved symmetries. 



KOFLiq=0 
KOFLiq=+l 
KOFLiq=-l 
K0FLiq=+2 
K0FLiq=-2 
K0FLiq=+3 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
KOFLiq 
K0FLiq=-9 
K0FLiq=+10 
K0FLiq=-10 



-3 
-4 
-4 
-5 
-5 
-6 
-6 
-7 
-7 



-9 



No diabatic blocking 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 
The state which has 



; in the given parity block. 

the larger ^/-alignment is occupied. 

the smaller ^/-alignment is occupied. 

the larger y-intrinsic spin is occupied. 

the smaller y-intrinsic spin is occupied. 

the larger x-alignment is occupied. 

the smaller x-alignment is occupied. 

the larger x-intrinsic spin is occupied. 

the smaller x-intrinsic spin is occupied. 

the larger z-alignment is occupied. 

the smaller ^-alignment is occupied. 

the larger z-intrinsic spin is occupied. 

the smaller ^-intrinsic spin is occupied. 

the larger o-alignment is occupied. 

the smaller o-alignment is occupied. 

the larger o-intrinsic spin is occupied. 

the smaller o-intrinsic spin is occupied. 

the larger multipole moment Q20 is occupied. 

the smaller multipole moment Q20 is occupied. 

the larger multipole moment Q22 is occupied. 

the smaller multipole moment Q22 is occupied. 



Here, the x-, y-, and ^-alignments or intrinsic spins denote projections of the total angular 
momentum or spin, respectively, on the x, y, and z axes. Similarly, the o-alignment or o-intrinsic 
spin denotes analogous projections on the direction of the angular frequency ujj. 

Within the diabatic blocking procedure one does not predefine whether the particle or the 
hole state is occupied (as is the case when the particle-hole excitations are defined, see Section 
3.4 of II). In each iteration the code calculates the average alignments (or average intrinsic 
spins, or average quadrupole moments) of both states (those defined by KPFLIM and KHFLIM), 
and occupies that state for which a larger, or a smaller value is obtained. Therefore, the order 
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of both states in the Routhian spectrum is irrelevant. 

The user is responsible for choosing the particle-state indices (in KPFLIM) only among those 
corresponding to empty single-particle states, and the hole-state indices (in KHFLIM) only among 
those corresponding to occupied single-particle states. 



2.6.6. Keyword: DIAPAR_PRO 

2,2, 1,1, 0,0= KPFLIQ(0,1), KPFLIQ(1,1), 

KHFLIQ(0,1), KHFLIQ(1,1), 
K0FLIQ(0,1), K0FLIQ(1,1) 

Same as above but for the diabatic blocking of proton single-particle parity configurations. 



2.6.7. Keyword: BLOCKSiq_N 

1, 0, = INSIQN, IPSIQN, IDSIQN 
For |IDSIQN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when parity is conserved. For IDSIQN=+1 or —1, the blocked quasiparticle state is selected 
by having the largest overlap with the INSIQN-th neutron single-particle eigenstate of the HFB 
mean-field Routhian in a given parity block or with its time-reversed partner, respectively. The 
parity of the state, +1 or —1, is defined by IPSIQN=0 or 1, respectively. |IDSIQN|=1 requires 
ISIMPY=0, IPARTY=1, IPAHFB=1, and IR0TAT=1. 



2.6.8. Keyword: BLOCKSIQ_P 

1, 0, = INSIQP, IPSIQP, IDSIQP 
Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 



2.7. Configurations with conserved parity and signature 
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2.7. Configurations with conserved parity and signature 

This section lists keywords pertaining to definitions of configurations in the case when parity 
and signature are conserved, and the single-particle or single-quasiparticle space is divided into 
eight charge-parity-signature blocks (neutrons and protons with parities 7r=±l and signatures 
r=±i). 



2.7.1. Keyword: VACSIG_NEU 

22, 22, 21, 21 =KVASIG(0,0), KVASIG(0,1), 
KVASIG(1,0), KVASIG(1,1) 

Numbers of lowest neutron states occupied in the four parity-signature blocks, denoted by 
(+,+), (+,—), (—,+), and (—,—), of given (paritysignature) combinations, i.e., (vr, r) = (+l,+i), 
(+1,-2), (— 1,+z), and (— 1,— i), respectively. These numbers define the parity-signature ref- 
erence configuration from which the particle-hole excitations are counted. The definitions of 
parity-signature reference configuration and excitations are ignored unless ISIMPY=1, ISIGNY= 
1, and IPAIRI=0. 



2.7.2. Keyword: VACSIG_PRO 



16, 16, 17, 17 =KVASIG(0,0), KVASIG(0,1), 
KVASIG(1,0), KVASIG(1,1) 



Same as above but for the numbers of proton states. 



2.7.3. Keyword: PHSIGN_NEU 

1, 00,00,00,00, 00,00,00,00 =NUPAHO, 

KPPPSP, KPPPSM, KPPMSP, KPPMSM, 
KHPPSP, KHPPSM, KHPMSP, KHPMSM 

Neutron particle-hole excitations in the parity-signature blocks. Basic principles are the same 
as for the excitations in the parity blocks. Particles are removed from the KHPPSP-th state in 
the (+,+) block, from the KHPPSM-th state in the (+,-) block, from the KHPMSP-th state in the 
(-,+) block, and from the KHPMSM-th state in the (-,-) block, and put in the KPPPSP-th state 
in the {+,+) block, in the KPPPSM-th state in the (+,-) block, in the KPPMSP-th state in the 
(— ,-|-) block, and in the KPPMSM-th state in the (— ,— ) block. These particle-hole excitations 
are ignored unless ISIMPY=1, ISIGNY=1, and IPAIRI=0. 



42 



2. INPUT DATA FILE 



2.7.4. Keyword: PHSIGN_PRO 

1, 00,00,00,00, 00,00,00,00 =NUPAHO, 

KPPPSP, KPPPSM, KPPMSP, KPPMSM, 
KHPPSP, KHPPSM, KHPMSP, KHPMSM 

Same as above but for the proton particle-hole excitations. 



2.7.5. Keyword: DIASIG_NEU 

2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, = 

KPFLIG (0,0,0), KPFLIG (0,1,0), KPFLIG(1 , , 0) , KPFLIG(1 , 1 ,0), 
KHFLIG (0,0,0), KHFLIG (0,1,0), KHFLIG(1,0,0), KHFLIG(1 , 1 ,0), 
KOFLIG(0,0,0), K0FLIG(0,1,0), K0FLIG(1 ,0,0), K0FLIG(1 , 1 ,0) 

The diabatic blocking of neutron single-particle parity-signature configurations. Matrices 
KPFLIG contain the indices of particle states in the four parity-signature blocks denoted by 
(-|-,+), (+,—), (—,+), and (—,—), of given (paritysignature) combinations, i.e., (vr,r) = 

and respectively. Matrices KHFLIG contain analogous indices of hole 

states, and matrices KOFLIG define the type of blocking according to Table HI Other rules de- 
scribed for the parity case apply here analogously. 



2.7.6. Keyword: DIASIG_PRO 

2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0, = 

KPFLIG (0,0,1), KPFLIG(0,1,1), KPFLIG(1,0, 1), KPFLIG(1 , 1 , 1), 
KHFLIG(0,0,1), KHFLIG(0,1,1), KHFLIG(1,0, 1), KHFLIG(1 , 1 , 1), 
K0FLIG(0,0,1), K0FLIG(0,1,1), K0FLIG(1 , , 1) , K0FLIG(1 , 1 , 1) 

Same as above but for the diabatic blocking of proton single-particle parity-signature configu- 
rations. 



2.7.7. Keyword: BLOCKSIG_N 

1, 0, 0, = INSIGN, IPSIGN, ISSIGN, IDSIGN 
For |IDSIGN|=1, the code performs the neutron quasiparticle blocking calculations in the case 
when parity and signature are conserved. For IDSIGN=-|-1 or —1, the blocked quasiparticle state 
is selected by having the largest overlap with the INSIGN-th neutron single-particle eigenstate of 
the HFB mean-field Routhian in a given parity-signature block or with its time-reversed part- 
ner, respectively. The parity of the state, +1 or —1, is defined by IPSIGN=0 or 1, respectively. 
The signature of the state, +i or —i, is defined by ISSIGN=0 or 1, respectively. |IDSIGN|=1 
requires ISIMPY=1, IPARTY=1, IPAHFB=1, and IR0TAT=1. 



2.7. Configurations with conserved parity and signature 
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2.7.8. Keyword: BLOCKSIG_P 

1, 0, 0, = INSIGP, IPSIGP, ISSIGP, IDSIGP 
Same as above but for the proton quasiparticle blocking. For odd-odd nuclei, neutron and 
proton quasiparticles can be simultaneously blocked. 
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2.8. Configurations 

This section lists keywords pertaining to definitions of configurations in the case when any 
combination of parity and signature symmetry is conserved. 

2.8.1. Keyword: BLOCKFIX_N 

0, = IFIBLN, INIBLN 

For IFIBLN=1, the neutron quasiparticle blocking is based on calculating overlaps with a fixed 
single-particle wave function. The method is based on beginning the iteration by defining the 
number of the single-particle state, that is, INSIGN, INSIMN, INSIQN, or INSIZN, depending on 
the selected symmetry conditions. Then, for INIBLN=1 this one single-particle wave function is 
stored in memory and on the record file, and in consecutive iterations the overlaps are calculated 
with respect to this fixed single-particle wave function. Therefore, subsequent changes in the 
ordering and structure of single-particle states do not affect the blocking mechanism. For 
INIBLN=0, in the first iteration this fixed single-particle wave function is not initialized, but 
it is read from the record file. IFIBLN=1 requires |IDSIGN|=1, |IDSIMN|=1, |lDSiqN|=l, or 
|IDSIZN|=1, depending on the selected symmetry conditions. IFIBLN=1 and INIBLN=0 requires 
IC0NTI=1. 



2.8.2. Keyword: BLOCKFIX_P 

0, = IFIBLP, INIBLP 
Same as above but for the proton quasiparticle blocking. 



2.9. Miscellaneous parameters 
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2.9. Miscellaneous parameters 

This section lists keywords pertaining to definitions of parameters specifying various numerical 
conditions of the code execution. 



2.9.1. Keyword: MAX_MULTIP 

2, 4, 4 = NMUCON, NMUCOU, NMUPRI 

Maximum multipolarities A of multipole moments used in the code for the constraints, Eq. f|T5|) . 
surface term of the Coulomb mean field, and printed on the output, respectively. Values not 
larger than A=9 are currently allowed. In case of the conserved parity, only even multipoles 
are used in the Coulomb mean field. 



2.9.2. Keyword: MAX_SURFAC 

0, = NSICON, NSIPRI 

Maximum multipolarities A of surface multipole moments, used in the code for the constraints, 
Eq. (fTSl) . and printed on the output, respectively. Values not larger than A=7 are currently 
allowed. 



2.9.3. Keyword: MAX_MAGNET 

0, = NMACON, NMAPRI 

Same as above but for the magnetic moments. In version (v2.40h) constraints on magnetic 
moments are not yet implemented, which requires NMACON=0. Values not larger than A=9 are 
currently allowed. 



2.9.4. Keyword: COULOMB 

80, 79, 0.25 = NUMCOU, NUMETA, FURMAX 

NUMCOU gives the number of points ]\['^°'^^, Eq. (I-lOO), used when summing up the Coulomb 
Green function. The dimensionless parameter d defining the size of the Coulomb parallelepiped, 
Eqs. (1-102) and (1-96), is given by NUMCOU* FURMAX. NUMETA gives the order of the Simpson 
integration of the solid harmonics on the faces of the parallelepiped, see Sec. 1-5. 



2.9.5. Keyword: SLOW-DOWN 

0.5, 0.5 = SLOWEV, SLOWOD 

The standard prescription to calculate the HP potential in the next iteration is to mix a given 
fraction e of the HP potentials from the previous iteration, with the fraction 1— e of potentials 
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calculated in the current iteration. SLOWEV and SLOWOD give the values of the mixing fractions e, 
separately for the time-even and time-odd potentials. This mixing slows down the convergence, 
and is necessary for the convergence of the iterative procedure. 



2.9.6. Keyword: SLOW-PAIR 

0.5 = SLOWPA 

Similarly as in the particle-hole channel, SLOWPA gives the value of the mixing fraction e used 
in the pairing channel, in analogy to SLOWEV and SLOWOD. 



2.9.7. Keyword: EPS_HERMIT 

l.OE-14 = EPSHER 

Numerical precision requested for determining values of Hermite polynomials. 



2.9.8. Keyword: OPTI_GAUSS 

1 = lOPTGS 

For I0PTGS=1 and IREAWS=0, expression = 2N^ + 2 is used to calculate the order of the 
Gauss-Hermite integration from the maximum number A'^^ of HO quanta in the direction 
/i=x, y, or z, and the input parameters NXHERM, NYHERM, and NZHERM are ignored. 



2.9.9. Keyword: GAUSHERMIT 

18, 18, 32 = NXHERM, NYHERM, NZHERM 

Orders L^., Ly, and of the Gauss-Hermite integration in three Cartesian directions. Must be 
even. Ignored if IOPTGS=0 or IREAWS=1. 



2.9.10. Keyword: BROYDEN 

0, 7, 0.8, 1000. = IBROYD, N_ITER, ALPHAM, BROTRI 
For IBR0YD=1, the Broyden method is used to accelerate the convergence, see Section VI-2.8. 
For IBROYD=0, or IBR0YD=1 and N_ITER=0, variables N_ITER, ALPHAM, and BROTRI, which are 
read on this line, are ignored, that is, the previous values are kept. N_ITER is the number of 
iterations used to approximate the inverse Jacobian, ALPHAM is the value of the parameter a of 
the linear mixing to which is added the Broyden correction, and BROTRI triggers an automatic 
switch to the Broyden method, that is, when the absolute value of the stability energy becomes 
lower than BROTRI, iterations are changed to the Broyden scheme. A large value of BROTRI 
ensures that the Broyden method is used from the very first iteration. IBR0YD=1 is incompatible 
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with I_YUKA>0, IC0UDI=2, or IC0UEX=2. The Broyden method is implemented only in the 
FORTRAN-90 version of hfodd. 
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2.10. Parameters of the HO basis 

This section lists keywords pertaining to definitions of parameters of the HO basis. 



2.10.1. Keyword: BASIS_SIZE 

15, 301, 800.0 = NOSCIL, NLIMIT, ENECUT 

The HO basis is composed of states having not more than A'"o=NOSCIL quanta in either of the 
Cartesian directions, and not more than M=NLIMIT states in totaL The states are added to 
the basis according to the increasing energy of the deformed harmonic oscillator, 

€-n^nyn, = ^^xin^ + \) + hUy{ny + |) + + \) . (12) 

In case of degenerate HO states (e.g., for an axially deformed HO) the complete multiplets 
are included, so the actual number of states LDBASE can be slightly larger than NLIMIT. If 
NLIMIT=0, it is set to the number of states up to NOSCIL spherical HO shells, that is, NLIMIT = 
((N0SCIL+l)*(N0SCIL+2)*(N0SCIL+3))/6. If NLIMIT>0, the cutoff energy ENECUT is ignored; 
otherwise all states having HO energy smaller then ENECUT are included in the basis, and NLIMIT 
is ignored. 



2.10.2. Keyword: HOMEGAZERQ 

1.2 = FCHOMO 

The code uses the standard value of the spherical HO frequency ^a;o=41 MeV/A^/^ multiplied 
by the scaling factor /=FCHOMO. 



2.10.3. Keyword: SURFAC_PAR 



86, 66, 1.23 



INNUMB, IZNUMB, ROPARM 



The code hfodd calculates parameters of the HO basis, and the zero-iteration Nilsson potential, 
by defining the standard nuclear surface S, 



E : /?(^,0) = c{a) 



1+ E E o^x,yx,{d,' 

A=0 /i=-A 



(13) 



corresponding to the volume |7r*R0PARM^*(lNNUMB+IZNUMB), where c{a) is a function of 
such that the volume enclosed by the surface S does not depend on a, i.e., in particular 
c(0) =i?o=ROPARM* ( INNUMB+ IZNUMB) . 



2.10. Parameters of the HO basis 
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2.10.4. Keyword: SURFAC_DEF 

2, 0, 0.61 = LAMBDA, MIU, ALPHAR 

The code defines frequencies of the deformed HO in three directions by using relation u^Rx = 
ojyRy = uJzRz, where i?^ are the lengths of principal axes of the nuclear surface S (fT3|) defined 
by real deformation parameters aA/x=ALPHAR(LAMBDA,MIU). The overall factor is defined by 

(uj^UJyUJzy/^=UJo. 
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2.11. Constraints on angular momentum 

This section lists keywords pertaining to definitions of constraints on angular momentum. The 
complete cranking term in the energy density functional is given by: 

/ c.,ox(Jo) V 



£ 



cran ^ ^ 



—^JOaiJoa) + Cja ({Joa) — Joa) — ^.Jla{Jla) 



2.11.1. Keyword: OMEGAY 

0.5 = OMEGAY 

Isoscalar angular frequency ujj^y in the y direction, see Eq. (fT4l) . A non-zero value of OMEGAY 
requires IR0TAT=1 and cannot be used for the BCS pairing, i.e., for IPAIRI=1 and IPAHFB=0. 



2.11.2. Keyword: OMISOY 

0.00 = OMISOY 

Same as above for the isovector angular frequency ujjiy in the y direction. 



2.11.3. Keyword: OMEGA_XYZ 

0.00, 0.00, 0.00, = OMEGAX, OMEGAY, OMEGAZ, ITILAX 

For ITILAX=1, values of the three Cartesian components of the isoscalar angular frequency 
vector cjjo, see Eq. (1141) . For ITILAX=0 these values are ignored. ITILAX=1 requires IR0TAT=1, 
ISIMPY=0, and IPAIRI=0. A non-zero value of OMEGAX, OMEGAY, or OMEGAZ requires broken 
symmetry ISIMTX=0, ISIMTY=0, or ISIMTZ=0, respectively. 



2.11.4. Keyword: OMEGA_RTP 

0.00, 0.00, 0.00, = OMERAD, OMETHE, OMEPHI, ITILAX 

Same as above but for values of the standard spherical components of the isoscalar angular 
frequency, i.e., wjor, c^joe, and lojq^. 



2.11.5. Keyword: OMISO_XYZ 

0.00, 0.00, 0.00, = OMISOX, OMISOY, OMISOZ, ITISAX 

For ITISAX=1, values of the three Cartesian components of the isovector angular frequency 
vector uJji, see Eq. ( !T4l) . For ITISAX=0 these values are ignored. ITISAX=1 requires IR0TAT=1, 
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ISIMPY=0, and IPAIRI=0. A non-zero value of OMISOX, OMISOY, or OMISOZ requires broken 
symmetry ISIMTX=0, ISIMTY=0, or ISIMTZ=0, respectively. 



2.11.6. Keyword: OMEGA_TURN 

= IMOVAX 

For IM0VAX=1 or —1, the isoscalar angular frequency vector ojjq is in each iteration set in the 
direction of the total angular momentum Jq, or opposite to this direction, respectively, while 
its length is kept equal to that fixed by the sum of squares of OMEGAX, OMEGAY, and OMEGAZ. For 
IMOVAX=0, vector ujjq is not changed during the iteration. IM0VAX=1 or —1 requires ITILAX=1, 
ITISAX=0, and IFLAGA=0. 



2.11.7. Keyword: SPINCONSTR 

0.0, 0.0, = STIFFI(2), ASKEDI(2), IFLAGI(2) 

For IFLAGI(2) = 1, the quadratic constraint on spin is used together with the linear constraint. 
Values of STIFFI(2) and ASKEDI(2) correspond respectively to Cjy and Joy in Eq. (fT4|) . For 
IFLAGI(2)=0, there is no quadratic constraint on spin. This keyword is maintained only 
for compatibility reasons; in version (v2.40h) of the code hfodd it has been replaced by 
SPICON_XYZ. 



2.11.8. Keyword: SPICON_XYZ 

0.0, 0.0, = STIFFI(l), ASKEDI(l), IFLAGI(l) 

0.0, 0.0, = STIFFI(2), ASKEDI(2), IFLAGI(2) 

0.0, 0.0, = STIFFKS), ASKEDKS), IFLAGI(3) 

Same as above but for quadratic constraints in the x, y, and z directions corresponding to 
IFLAGI(1)=1, IFLAGI(2)=1, and IFLAGI(1)=1, respectively. IFLAGI(1)=1, IFLAGI(2)=1, 
or IFLAGI(3)=1 requires broken symmetry ISIMTX=0, ISIMTY=0, or ISIMTZ=0, respectively. 



2.11.9. Keyword: SPICON_OME 

0.0, 0.0, = STIFFA, ASKEDA, IFLAGA 

For IFLAGA=1, the quadratic constraint on the angle between the angular frequency and angular 
momentum vectors, see Eq. ( IT^ . In version (v2.40h) the angle is constrained to zero. Value 
of STIFFA corresponds to Ca in Eq. f|T^ . Value of ASKEDA must be set to 0; this variable is 
reserved for a future implementation of the constraint to a non-zero angle. For IFLAGA=0, there 
is no quadratic constraint on the angle. IFLAGA=1 requires ISIMPY=0, IR0TAT=1, ITILAX=1, 
ITISAX=0, IFLAGI(1)=0, IFLAGI(2)=0, IFLAGI(3)=0, and a non-zero value of sum of squares 
of OMEGAX, OMEGAY, and OMEGAZ. 
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2.11.10. Keyword: NORBCONSTR 

= NO_ORB 

For N0_0RB=1, constraints on the intrinsic spin only are used and the orbital part of the angular 
momentum is not constrained. N0_0RB=1 requires IR0TAT=1. 



2.12. Constraints on multipole moments 
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2.12. Constraints on multipole moments 

This section lists keywords pertaining to definitions of constraints on multipole moments and 
surface multipole moments. The complete multipole constraint term in the energy density 
functional is given by: 

The multipole and surface multipole moments are defined as 

QA/i(r)=a,^rX(^>0). (16) 

and 

QUr) = a,,r'^'Y*^{e,<j,), (17) 

respectively, where Yxf^, are the standard spherical harmonics in the convention of Ref . [ID] • For 
A<2, factors are listed in Tabled while for A>2 axij,=l- 

Table 5: Adopted definitions of the normalization factors ax^ and the corresponding multipole 
moments Qxfi- 



X ji 


O-A/i 


Qa/1 







1 


1 


\/4vr/3 


z 


1 1 




X — iy 


2 


yi67r/5 


2^2 _ a;2 _ ^2 


2 1 


-^Stt/IS 


zx — izy 


2 2 


y327r/5 


V3(x^ — — 2ixy) 



2.12.1. Keyword: MULTCONSTR 

2, 0, 0.01, 42.0, 1 = LAMBDA, MIU, STIFFQ, QASKED, IFLAGQ 

For IFLAGQ=1, the total multipole moment of the given multipolarity A and /i is constrained. 
Values of LAMBDA, MIU, STIFFQ, and QASKED correspond respectively to A, /i, Ca/j, and Qa/i in 
Eq. (fT5|) . For IFLAGQ=0, there is no constraint in the given multipolarity. 



2.12.2. Keyword: SURFCONSTR 

2, 0, 0.0, 0.0, = LAMBDA, MIU, STIFFS, SASKED, IFLAGS 
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For IFLAGS=1, the total surface multipole moment of the given muhipolarity A and /i is con- 
strained. Values of LAMBDA, MIU, STIFFS, and SASKED correspond respectively to A, /i, Cf^, and 
Qf^ in Eq. For IFLAGS=0, there is no constraint in the given multipolarity. For IFLAGS=2 

or 3, only the neutron or proton, respectively, surface multipole moments are constrained. For 
all constrained multipolarities, values of IFLAGS must be the same. 



2.12.3. Keyword: MULTCONSCA 

1, 0., 0., = LAMBDA, STIFFG, GASKED, IFLAGG 

For IFLAGG=1, the total scalar multipole moment of the given multipolarity A is constrained. 
Values of LAMBDA, STIFFG, and GASKED correspond, respectively, to A, Cx, and Qx in Eq. (VI- 
80). For IFLAGG=0, there is no scalar constraint in the given multipolarity. The constrained 
multipolarity LAMBDA cannot be equal to or exceed NMUCON. For conserved parity IPARTY=1, 
only even moments can be constrained. 



2.12.4. Keyword: BQHR_BETAS 

4, 0, 1 = NEXBET, IPRIBE, IPRIBL 
For a given set of electric multipole moments, the code calculates and prints the corresponding 
first-order (for IPRIBL=1) and/or exact (for IPRIBE=1) Bohr deformation parameters, see Sec- 
tion VI-2.5. For IPRIBE=1, in case the code fails to find the exact values, the first-order values 
are printed irrespective of the value of IPRIBL. The exact values are sought for multipolarities 
up to NEXBET, which must not be greater than NMUPRI. The approximate values are printed up 
to multipolarity of NMUPRI. 



2.12.5. Keyword: SCHIFF_MOM 

= ISCHIF 

For ISCHIF=1, the surface multipole moments are everywhere in the code replaced by the Schiff 
moments, see Section VI-2.6. 



2.13. Output-file parameters 
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2.13. Output-file parameters 

This section lists keywords pertaining to definitions of parameters specifying character of in- 
formation printed on the output file. 



2.13.1. Keyword: PRINT-ITER 

1, 0, 1 = IPRSTA, IPRMID, IPRSTO 

The code prints results for the first, middle, and/or last iteration if the corresponding param- 
eters IPRSTA, IPRMID, and IPRSTO equal 1. 



2.13.2. Keyword: PRINT-MOME 

1, 1, 1 = IPRI_N, IPRI_P, IPRI_T 

The code prints values of the neutron, proton, and total moments (multipole moments, surface 
multipole moments, or magnetic moments) only if the corresponding parameters IPRI_N=1, 
IPRI_P=1, and IPRI_T=1, respectively. 



2.13.3. Keyword: PRINT-INTR 

1 = INTRIP 

The code prints the values of moments (multipole moments, surface multipole moments, or 
magnetic moments) and angular momenta in the intrinsic frame only if INTRIP=1. 



2.13.4. Keyword: EALLMINMAX 

-12.0, 0.0 = EMINAL, EMAXAL 

The code prints tables of single-particle properties for states with values of the single-particle 
Routhians between EMINAL and EMAXAL. No table is printed unless EMINAL < EMAXAL. 



2.13.5. Keyword: EQUASI_MAX 

10.0 = EMAXQU 

The code prints tables of quasiparticle properties for the HFB states with values of the quasi- 
particle energies smaller than EMAXQU. No table is printed unless EMAXQU>0. 
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2.13.6. Keyword: HFBMEANFLD 

= IMFHFB 

For IMFHFB=1, eigenvalues of the HFB mean- field single-particle Hamiltonian or Routhian are 
printed. IMFHFB=1 requires IPAHFB=1. 



2.13.7. Keyword: PRINT_AMP 

0, 999, 1, 1, 1, 0, 

ISLPRI, ISUPRI, lENPRI, ISRPRI, IMIPRI, IKEPRI, IRMPRI 

These parameters govern the printing of the AMP results. If not restricted by the values of 
the doubled angular momenta for which the calculations are performed, IPROMI and IPROMA, 
the results are printed only for the doubled angular momenta between ISLPRI and ISUPRI. For 
IENPRI=1, the AMP energies are printed. In addition, for IENPRI=2, the AMP kernels are also 
printed. For ISRPRI=1, the sum rules are printed and compared with the HF average values. 
For IMIPRI=1, the energies of the i^-mixed states are printed. For IKEPRI=1 and/or IRMPRI=1, 
the proton reduced kernels and/or reduced matrix elements, respectively, are printed. 



2.13.8. Keyword: TRANCUTPRI 

0., 0., 0., = QMUCUT, QMACUT, QSICUT 
Values of the proton kernels and reduced matrix elements for the transition electric, transition 
magnetic, and transition surface or Schiff moments are printed only if their absolute values are 
larger than, respectively, QMUCUT, QMACUT, and QSICUT. This option is meant to avoid printing 
long lists of very small or zero values. 



2.13.9. Keyword: ONE_LINE 

1 = IlLINE 

For IILINEt^ 0, a one-line convergence report is printed at each iteration. For I1LINE=1, the 
code prints the values of deformation 7, total angular momentum, total angular frequency u;, 
and angle between vectors of angular momentum and frequency. For I1LINE=2, the code prints 
the values of neutron and proton pairing gaps and Lipkin-Nogami parameters A2. 



2.13.10. Keyword: PRINT_VIOL 

= IVIPRI 

For IVIPRI=1, the code prints integrals of several symmetry-violating terms. IVIPRI=1 re- 
quires IR0TAT=1. 



2.13. Output-file parameters 
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2.13.11. Keyword: NILSSONLAB 

3 = NILXYZ 

This option used for NILXYZ=1, 2, or 3, allows for printing the Nilsson labels defined with 
respect to the x, y, or z axis, respectively. This feature is useful when the axial symmetry axis 
of a nucleus and its angular momentum are aligned along the x ot y axis and not along the z axis, 
for which the standard Nilsson labels are defined. This is particularly important for analyzing 
configurations of band heads of odd nuclei within the conserved ^/-signature-symmetry limit. 
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2.14. Files 

This section lists keywords that define files used by the code. 
2.14.1. Keyword: REVIEWFILE 

HFODD.REV = FILREV 

CHARACTER*68 file name of the review file. Must start at the 13-th column of the data line. 
The ASCII review file is written after calculating every data set (i.e., once per every EXECUTE 
item), provided IREVIE>0. The file is not rewound, so the results for several data sets can be 
accumulated in a single file. This is so provided the filename FILREV is not changed between 
the EXECUTE items. The file contains all relevant parameters and results of calculation in a 
form suitable for reading by another program. It is meant as an interface to programs which 
analyze and/or plot the results. The file contains sections defined by keywords (different than 
keywords used in the input data file described here). The detailed structure of the review file 
is not documented in the present write up, and can be inferred from inspecting the specimen 
produced by the sample run, and from the subroutine review. 



2.14.2. Keyword: REVIEW 

2 = IREVIE 

The review file is not written if IREVIE=0. For IREVIE=2, the table of single-particle properties 
is included in the review file in addition to other results, which are written for IREVIE=1. 



2.14.3. Keyword: REPLAYFILE 

HFODD.REP = FILREP 

CHARACTER*68 file name of the replay file. Must start at the 13-th column of the data line. The 
binary replay file with the name defined in FILREP must exist if IC0NTI=1, and will be read. If 
the filenames FILREP and FILREC are identical, the replay file will be subsequently overwritten 
as a new record file. These feature is implemented to facilitate chaining of jobs which follow 
one another. 



2.14.4. Keyword: RECORDFILE 

HFODD.REC = FILREC 

CHARACTER*68 file name of the record file. Must start at the 13-th column of the data line. If 
IWRIRE=1, binary record file is written after each HF iteration. It contains complete information 
which allows restarting the iteration in another run of the code. To restart, one has to specify 
IC0NTI=1 and provide the name of the file by defining FILREP. In case of the computer crash, 
or upon a successful completion of the given input data set, the record file contains the results 
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of the last performed iteration. The file is always rewound before it is written, so the results 
for consecutive iterations do not pile up. 



2.14.5. Keyword: RECORDSAVE 

1 = IWRIRE 

For IWRIRE=1, the record file is saved on disc after each iteration is completed. For IWRIRE=0, 
it is saved only once, after all iterations are completed. This option is useful on systems where 
the speed of writing large amounts of data to disc could hamper the system performance. For 
IWRIRE=— 1, the record file is never saved, and the run cannot be later continued. 



2.14.6. Keyword: WOODSAFILE 

WOODS.WFN = FILWOO 

CHARACTER*68 file name of the Woods-Saxon file. Must start at the 13-th column of the 
data line. This file is read provided IREAWS=1. The binary Woods-Saxon file constitutes an 
interface between the Woods-Saxon code and the code hfodd. It contains the Woods-Saxon 
wave functions and numerous other parameters which define the current calculation. The 
parameters read from the Woods-Saxon file overwrite the values provided in the input data file. 
Since the current version of the Woods-Saxon code will be published separately, the feature of 
starting the iteration from the Woods-Saxon results is not documented in the present write up. 



2.14.7. Keyword: COULOMFILE 

HFODD. COU = FILCOU 

CHARACTER*68 file name of the Coulomb file containing auxiliary data for a faster calculation of 
the Coulomb mean field. Must start at the 13-th column of the data line. Parameters NUMCOU, 
NUMETA, and FURMAX which define the calculations of the Coulomb direct potential (see Sees. 
1-5 and II-3.5) are usually kept unchanged for the whole series of calculations in one region of 
nuclei. Therefore, many Coulomb auxiliary results can be calculated only once, and stored in 
the file FILCOU. Handling of this file is determined by the input parameters ICOULI and ICOULO. 



2.14.8. Keyword: COULOMSAVE 

0,0 = ICOULI, ICOULO 

Input parameters ICOULI and ICOULO determine actions pertaining to reading and/or writing 
of the Coulomb file FILCOU, according to the following table: 
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ICOULI ICOULO Action 

neither read nor write the Coulomb file 

1 read, but do not write the Coulomb file 

1 do not read, but write the Coulomb file 

1 1 use automated handling of the Coulomb file 



The default values of ICOULI=ICOULO=0 ensure that the code hfodd (v2.40h) behaves as 
that in version (vl.60r). However, unless it is required by special circumstances, a use of the 
automated handling of the Coulomb file (IC0ULI=IC0UL0=1) is recommended. Within the 
automated mode, the code checks whether the Coulomb file exists, and whether it contains 
data which match the current values of the input parameters NUMCOU, NUMETA, and FURMAX. If 
this is the case, the code reads the data from the Coulomb file. If this is not the case, the 
code calculates the Coulomb auxiliary results and stores them in the Coulomb file FILCOU. In 
the automated mode, the user is informed by appropriate messages printed on the output file 
about what type of the action has been taken in the given run of the code. 



2.14.9. Keyword: WAVEF_FILE 

HFODD.WFN = FILWAV 
CHARACTER*68 file name of the wave function file. Must start at the 13-th column of the data 
line. The binary wave function files must exist if IPRGCM=2, and will be read, see the keyword 
PROJECTGCM. 



2.14.10. Keyword: KERNELFILE 

HFODD. KER = FILKER 
CHARACTER*68 file name of the kernel file. Must start at the 13-th column of the data line. The 
binary kernel files are written and read if ISAKER=1, see the keyword SAVEKERNEL. 



2.14.11. Keyword: YUKAWASAVE 

1 = IWRIYU 

For IWRIYU=1 and I_YUKA>2, the Yukawa file is saved on disc after each iteration is completed. 
The file contains the matrix elements of the Yukawa mean field. For IWRIYU=0 and I_YUKA>2, 
the file is saved only once, after all iterations are completed. For IWRIYU=— 1, the file is 
never saved. The Yukawa calculations can be restarted by using the Yukawa file. Therefore, 
IYC0NT=1 requires setting IWRIYU=0 or 1 in the run that is to be continued. Equivalently, the 
Yukawa calculations can be restarted by using the field file, which requires setting IWRIFI=0 
or 1 in the run that is to be continued. 



2.14.12. Keyword: REPYUKFILE 

HFODD.YUP = FILYUP 
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CHARACTER*68 file name of the Yukawa file. Must start at the 13-th column of the data line. 
The binary Yukawa file with the name defined in FILYUP must exist if IYC0NT=1, and will be 
read. If the filenames FILYUP and FILYUC are identical, the Yukawa file will be subsequently 
overwritten as a new Yukawa file. 



2.14.13. Keyword: RECYUKFILE 

HFODD.YUC = FILYUC 

CHARACTER*68 file name of the Yukawa file. Must start at the 13-th column of the data line. 
If IWRIYU=1, binary Yukawa file is written after each HP iteration. It contains complete 
information that allows restarting the Yukawa calculations in another run of the code. To 
restart, one has to specify IYC0NT=1 and provide the name of the file by defining FILYUP. 



2.14.14. Keyword: LIPKINSAVE 

1 = IWRILI 

For IWRILI=1 and LIPKIN=1 or LIPKIP=1, the Lipkin-Nogami file is saved on disc after each 
iteration is completed. The file contains the matrix elements of the particle density matrices. 
For IWRILI=0 and LIPKIN=1 or LIPKIP=1, the file is saved only once, after all iterations are 
completed. For IWRILI=— 1, the file is never saved. The Lipkin-Nogami calculations can be 
restarted by using the Lipkin-Nog ami file. Therefore, ILC0NT=1 requires setting IWRILI=0 
or 1 in the run that is to be continued. Equivalently, the Lipkin-Nogami calculations can be 
restarted by using the field file, which requires setting IWRIFI=0 or 1 in the run that is to be 
continued. 



2.14.15. Keyword: REPLIPFILE 

HFODD.LIP = FILLIP 

CHARACTER*68 file name of the Lipkin-Nogami file. Must start at the 13-th column of the data 
line. The binary Lipkin-Nogami file with the name defined in FILLIP must exist if ILC0NT=1, 
and will be read. If the filenames FILLIP and FILLIC are identical, the Lipkin-Nogami file will 
be subsequently overwritten as a new Lipkin-Nogami file. 



2.14.16. Keyword: RECLIPFILE 

HFODD.LIC = FILLIC 

CHARACTER*68 file name of the Lipkin-Nogami file. Must start at the 13-th column of the data 
line. If IWRILI=1, binary Lipkin-Nogami file is written after each HF iteration. It contains 
complete information that allows restarting the Lipkin-Nogami calculations in another run of 
the code. To restart, one has to specify ILC0NT=1 and provide the name of the file by defining 
FILLIP. 
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2.14.17. Keyword: FIELD_SAVE 

-1 = IWRIFI 

For IWRIFI=1, the field file, is saved on disc after each iteration is completed. The file contains 
the matrix elements of the mean field. For IWRIFI=0, the file is saved only once, after all 
iterations are completed. For IWRIFI=— 1, the file is never saved. To restart calculations of the 
exact Coulomb exchange energy, the field file must exist. IFC0NT=1 requires setting IWRIFI=0 
or 1 in the run that is to be continued. 



2.14.18. Keyword: FIELD_OLD 

= IWRIOL 

For IWRI0L=1, in the last iteration the mean fields are not updated, that is, the mixing fractions 
e are set equal to 1, irrespective of values of the SLOWEV, SLOWOD and SLOWPA parameters. In this 
way, information stored on the record file corresponds to the last but one iteration. Only then, 
restarting of the calculation from the field file leads to a smooth continuation of iterations. 



2.14.19. Keyword: REP_FIELDS 

HFODD.FIP = FILFIP 

CHARACTER*68 file name of the field file. Must start at the 13-th column of the data line. The 
binary field file with the name defined in FILFIP must exist if IFC0NT=1, and will be read. If 
the filenames FILFIP and FILFIC are identical, the field file will be subsequently overwritten 
as a new field file. 



2.14.20. Keyword: REC_FIELDS 

HFODD.FIC = FILFIC 

CHARACTER*68 file name of the field file. Must start at the 13-th column of the data line. If 
IWRIFI=1, binary field file is written after each HF iteration. It contains complete information 
that allows restarting the calculations from the matrix elements of fields. To restart, one has 
to specify IFC0NT=1 and provide the name of the file by defining FILFIP. 



2.15. Starting the iteration 
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2.15. Starting the iteration 

This section lists keywords that pertain to starting and restarting the code. 



2.15.1. Keyword: RESTART 

= ICONTI 

For IC0NTI=1, results stored in the replay file (written as a record file in a previous run) are 
used to start the iteration. The code assumes (without checking) that the parameters of the 
HO basis are the same as those used in the previous run. The replay file name should be 
provided by defining FILREP. If the previous run was done with IREAWS=1, the current run 
must also use IREAWS=1, and the same Woods-Saxon file must be provided. This so because 
the Woods-Saxon file contains not only the information about the starting potential (which is 
ignored for IC0NTI=1) but also defines the HO basis. 



2.15.2. Keyword: CONT_PAIRI 

= IPCONT 

For IPC0NT=1, results stored in the replay file are used to define the HFB pairing properties in 
the first iteration; otherwise values read in matrices FERINI and DELINI are used. When the 
HFB pairing correlations are taken into account (i.e., for IPAHFB=1), and if a smooth restart 
and continuation of iterations from previously stored results is required, value of IPC0NT=1 
must be used. IPC0NT=1 is incompatible with either of IPAHFB=0 or ICONTI=0. 



2.15.3. Keyword: CONT_OMEGA 

= lOCONT 

For I0C0NT=1, the angular momentum vector Jq stored in the replay file is used to define in the 
first iteration the direction of the angular frequency vector uj jo, see Eq. f|T^ : otherwise values 
read in OMEGAX, OMEGAY, and OMEGAZ are used (see keyword OMEGA_XYZ). I0C0NT=1 requires 
IM0VAX=1 or -1. 



2.15.4. Keyword: READ_WOODS 

= IREAWS 

For IREAWS=1, the results stored by the previously performed Woods-Saxon calculation are 
used to start the iteration and to define the HO basis. For IC0NTI=1 and IREAWS=1 the 
Woods-Saxon file must also be provided, and is used only to define the HO basis. 
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2.15.5. Keyword: NILSSONPAR 

0, -1.175, -0.247, -1.175, -0.352, 11.17, 11.17, 6.28 = 
NILDAT, CNILSN, DNILSN, CNILSP, DNILSP, 
HBANIX, HBANIY, HBANIZ 

For IREAWS=0 and ICONTI=0, the code starts the calculation from the Nilsson potential. For 
NILDAT=1 or NILDAT=2, the Nilsson parameters C and D [Ref. [II], Eq. (2.89)] are for neutrons 
given by CNILSN and DNILSN and for protons by CNILSP and DNILSP. For NILDAT=1, frequencies 
of the deformed HO are defined by HBANIX, HBANIY, and HBANIZ, while for NILDAT=2 they are 
the same as those of the HO basis defined under keyword SURFAC_DEF. If NILDAT=0, the Nilsson 
parameters C and D are defined by Eq. (2.91) and Table 2.3 of Ref. [H], while frequencies of 
the deformed HO are the same as those of the HO basis defined under keyword SURFAC_DEF. 
For NILDAT=0, the Nilsson parameters read from the input file are ignored. 



2.15.6. Keyword: CONTYUKAWA 

= lYCONT 

For IYC0NT=1, results stored in the Yukawa file are used to define the Yukawa fields in the 
first iteration; otherwise the Yukawa fields are set equal to zero. When the Yukawa fields are 
taken into account (i.e., for I_YUKA=2), and if a smooth restart and continuation of iterations 
from previously stored results is required, value of IYC0NT=1 must be used. IYC0NT=1 is 
incompatible with either of I_YUKA<2 or ICONTI=0. 



2.15.7. Keyword: CONTLIPKIN 

= ILCONT 

For ILC0NT=1, results stored in the Lipkin-Nogami file are used to define the density matrix 
required for the Lipkin-Nogami calculations in the first iteration; otherwise the required den- 
sity matrix is set equal to zero. For the Lipkin-Nogami calculations (i.e., for LIPKIN=1 or 
LIPKIP=1), and if a smooth restart and continuation of iterations from previously stored re- 
sults is required, value of ILC0NT=1 must be used. ILC0NT=1 is incompatible with either of 
LIPKIN=LIPKIP=0 or IPCONT=0. 



2.15.8. Keyword: CONTFIELDS 

= IFCONT 

For IFC0NT=1, results stored in the field file are used to define the matrix elements of fields in 
the first iteration; otherwise the matrix elements are recalculated. For the Gaussian-expansion 
method used to calculate the Coulomb energy and Coulomb mean field, that is for IC0UDI=2 or 
IC0UEX=2, and if a smooth restart and continuation of iterations from previously stored results 
is required, value of IFC0NT=1 must be used. IFC0NT=1 is incompatible with ICONTI=0. 



3. OUTPUT FILE 
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3. OUTPUT FILE 

The code writes results on standard FORTRAN output file. The output file begins with the 
information pertaining to the general parameters of the calculation, then gives information 
about the starting point of the iteration, provides the convergence report, and finally contains 
the results calculated at the last iteration. Structure of the output file strongly depends on the 
type of calculation that are requested in the given run; especially on the assumed conserved 
and broken symmetries. As a rule, only the information relevant for the specific run is printed. 
Descriptions and comments printed on the output file are in most cases self-explanatory. In 
this section we list only those explanations that are not explicitly included in the output file. 

• Section CLASSICAL NUCLEAR SURFACE lists deformation parameters used to define the 
nuclear surface from which the basis parameters are derived, see keyword SURFAC_DEF. It 
also gives the basis parameters such as the OSCILLATOR FREQUENCIES : HBAROX, HBAROY, 
and HBAROZ corresponding to TiuOx-, huy, and huz- 

• Section PHYSICAL CONSTANTS gives values of he in MeVfm (H_BARC), of hc/e"^ (HBC0E2), 
of the neutron and proton masses in MeV/c^ (XMASSN and XMASSP), of the kinetic -energy 
coefficient /2m before (HBMASS) and after the center-of-mass correction (HBMRPA), of 
the elementary charge squared in MeVfm (ECHAR2), and of the coefficient preceding 
the integral in the Coulomb exchange energy (COULEX). For details see the comments in 
the SETBAS subroutine. 

• Section OSCILLATOR LENGTHS, CONSTANTS, and FREQUENCIES gives the values of 1/6^, 6^, 
and huj^, respectively, which characterize the HO basis in three Cartesian directions, ^=x, 

• Section BASIS CUT-OFF CONTROL PARAMETERS gives maximum numbers of the HO quanta 
in three directions NXMAXX, NYMAXX, and NZMAXX corresponding to N^, Ny, and A''^, as well 
as the orders of the Gauss- Hermite quadratures NXHERM, NYHERM, and NZHERM, correspond- 
ing to La;, Ly, and L^, see keyword OPTI_GAUSS. It also gives the number M of the HO 
states included in the basis as requested in the input file (NLIMIT) and as used in the 
calculation (LDBASE). 

• Section SHAPE OF THE OSCILLATOR-BASIS DIAMOND gives the numbers of the HO quanta 
in a given direction for fixed numbers of the HO quanta in both remaining directions. 
The output is arranged in such a way that the shape of the grid of points n^nyUz is clearly 
visualized by projections in every of the three directions. 

• Section SKYRME FORCE DEFINITION lists the name and parameters of the Skyrme force, 
together with parameters KETA_J, KETA_W, KETACM, and KETA_M that define the way the 
given force should be used (see keyword SKYRME-STD). 

• Section COEFFICIENTS DEFINING THE SKYRME FUNCTIONAL gives values of the coupling 
constants in the Skyrme functional. These values take into account scaling factors, which 
are printed in the section SCALING FACTORS, unless they are all equal to 1. 
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• Section CALCULATIONS WITH THE TILTED-AXIS CRANKING gives values of components of 
the angular frequency vector and its length. For switches IM0VAX=1 or —1 and IOCONT=0, 
these values are used only in the first iteration, and later ignored, because the angular 
frequency vector is in each iteration readjusted to be aligned or anti-aligned with the 
angular momentum vector. 

• Sections PARITY/SIGNATURE, SIMPLEX, or PARITY CONFIGURATIONS give the vacuum and 
particle-hole configurations requested in the input file for the corresponding symmetry. 

• Section CONVERGENCE REPORT gives the list of performed iterations. For each iteration, 
one line is printed showing the energy, stability energy, average values of the intrinsic 
quadrupole moment Q and 7 deformation, total angular momentum, angular frequency, 
angle between the angular frequency and angular momentum vectors, and total pairing 
energy. 

• Section SINGLE-PARTICLE PROPERTIES lists the single-particle states calculated for the 
Nilsson, Woods-Saxon, or Hartree-Fock Routhian operators. For every state one line is 
printed which gives the value of the single-particle Routhian, the consecutive numbers 
in the parity-signature block or simplex block, the quantum numbers [A^, n2,A]f2 of the 
Nilsson asymptotic state which has the largest component in the given state, the average 
value of the parity operator (in %), the average values of projections of the intrinsic 
and total angular momenta (in Ti) , and their ratio called the (7- factor. For broken simplex 
symmetry, three projections of the total angular momentum and intrinsic spin are printed 
in the first and second line, respectively, for each single-particle state. 

• Section DENSITY INTEGRALS IN THE SKYRME FUNCTIONAL gives integrals of products of 
densities, which appear in the Skyrme functional. Terms in the functional are identified 
by the acronyms described under keywords EVE_SCA_TS and ODD_SCA_TS. 

• Section CONTRIBUTIONS TO ENERGY IN THE SKYRME FUNCTIONAL gives the values of var- 
ious terms which appear in the Skyrme functional. These contributions are the products 
of the coupling constants and of the density integrals described above. The sums of 
time-even and time-odd contributions are also printed. 

• Section EULER ANGLES OF THE PRINCIPAL-AXES FRAME gives the standard Euler angles 
a, /3, and 7 [lO] in degrees, which define the orientation of the intrinsic frame (principal- 
axes frame) of reference with respect to the original frame. 

• Section MULTIPOLE MOMENTS gives values of neutron, proton, or mass multipole moments 
with respect to the original frame of reference. In particular, the A=0 moment corre- 
sponds to the number of particles. Similarly, section MULTIPOLE MOMENTS [UNITS: (10 
FERMI) "LAMBDA] [INTRINSIC FRAME] gives analogous values with respect to the intrinsic 
frame of reference. Whenever the parity symmetry is broken, yet another analogous sec- 
tion gives information on values of multipole moments with respect to the center-of-mass 
reference frame. 

• Sections ROOT-MEAN-SQUARE AND GEOMETRIC SIZES gives the rms average values of the 
radius and of the x, and z coordinates. In order to better visualize the size of the 
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nucleus, the geometric sizes are also calculated by multiplying the rms radius by yS/S 
and the rms coordinates by y^- 

• Sections MAGNETIC MOMENTS [MAGNETON*FERMI" (LAMBDA-1)] and MAGNETIC MOMENTS 
[MAGNET0N*FERMI"(LAMBDA-1)] [INTRINSIC FRAME] give values of magnetic moments 
with respect to the original and intrinsic frame, respectively. 

• Section ANGULAR MOMENTA gives for ISIMPY=1 the average values of the total and intrin- 
sic neutron, proton, and total angular momentum (in h). It also gives the correspond- 
ing values and contributions to the first moment of inertia j'^^^=I /uj (in /i^/MeV). For 
ISIMPY=0, this section gives values of neutron, proton, and total projections of the total 
angular momentum and intrinsic spin on the three Cartesian axes in the original and 
intrinsic frame, respectively. 

• Sections NEUTRON CONFIGURATIONS and PROTON CONFIGURATIONS give a visual represen- 
tation of states occupied in the parity-signature blocks, simplex blocks, or parity blocks. 
The lines denoted by CONF: give the configurations requested in the input data set, while 
those denoted by VACC: give the configurations characterizing the given HF state. By 
comparing the two sets one can verify whether the requested configuration has been 
obtained, and eventually devise a new configuration to be calculated. The consecutive 
numbers printed in the horizontal direction correspond to the consecutive numbers in 
blocks printed in the section SINGLE-PARTICLE PROPERTIES. By comparing the two sec- 
tions one can effectively associate the Nilsson labels with the calculated configurations 
and also prepare the configuration input data. 

• Section ENERGIES gives a summary of the energies calculated for the HF state. The kinetic 
energy, single-particle energy, and pairing energy are printed for neutrons (NEU) , protons 
(PRO), and all particles (TOT). Then, pairing rearrangement energy P-REARR, pairing 
gap PAIRGAP, and Fermi energy E-FERMI are printed. The Coulomb direct energy (DIR) 
and Coulomb exchange energy (EXC) are printed together with their sum (TOT). The 
multipole (MULT), surface multipole (SURF), and cranking (SPIN) constraint energies, 
are printed together with the corresponding corrections energies (CORR.). Then the 
rearrangement energy is printed followed by the values of the Routhian and the spin-orbit 
and Skyrme energies, the latter two split in the time-even (EVE) and time-odd (ODD) 
contributions. Finally, the total energies are printed as (SP) and (FUN), respectively, 
while their difference is printed as the stability (STAB). 

• Section CODE COMPILED WITH THE... lists the values of the most important array dimen- 
sions and switches declared in the PARAMETER statements, fixed at the compilation stage. 

• Sections BOHR DEFORMATIONS (EXACT MULTIPOLE MOMENTS) give values of the Bohr de- 
formation parameters, determined as described in Section VI-2.5. 

• Section KERNELS AND AVERAGE VALUES... gives values of the diagonal norm kernels iV^i^ 
(VI-5) and average AMP energies Hkk/Nkk for different values of / and K. Real and 
imaginary parts of Nkk are both printed, although all imaginary parts should be equal 
to zero when the numerical precision is sufficiently good. 
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• Section SUM RULES... compares real and imaginary parts of the sum rules (VI-49) and 
(VI-50) with the corresponding HF average values. Results are printed for the norm 
kernels N^p^, which must sum up to 1. This condition is a primary test of whether 
a sufficient number of the angular momenta are included in the sum of Eq. (VI-48). 
In the example of the file ge064-a.out, the sum rule for the norm kernel NORM= 0.3248 
indicates a much too small value of Imax and too small numbers of the Gauss-Tchebyschev 
and Gauss-Legendre integration nodes. 

• Sections REDUCED KERNELS... give values of the reduced kernels (VI-56) of electric and 
magnetic transition operators. 

• Sections RESULTS OF THE K-MIXING... give values of the energies Ei of A'-mixed states, 
see Eq. (VI-3) and values of norm eigenvalues n^, see Eq. (VI-9). 

• Sections REDUCED MAT.ELEMS... give values of the reduced matrix elements (VI-57) of 
electric and magnetic transition operators calculated for the i^-mixed states. 

• Section NUMBERS OF CALLS TO SUBROUTINES gives the statistics of calls to subroutines, 
which together with the section EXECUTION TIMES IN SUBROUTINES illustrates the work 
flow of the code. 
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FORTRAN source of code hfodd (v2.40h) is provided in the file hfodd.f and can be modified 
in several places, as described in this section. 

4.1. Dimensions of arrays 

The code hfodd uses the arrays' dimensions declared through the PARAMETER statements. This 
allows changing the dimensions and adapting the size of the reserved memory to the problem 
being solved. Whenever too small a dimension is defined, the code aborts with a message 
indicating the dimension that should be increased. Substantial amount of memory is required 
only for arrays depending on the following PARAMETER values: 

PARAMETER (NDMAIN=21) 

Corresponds to the maximum number of the HO quanta. Should be larger or equal to the input 
parameter NOSCIL defined under keyword BASIS_SIZE. 

PARAMETER (NDBASE=507) 

Corresponds to the maximum number of the HO basis states. Should be larger or equal to the 
input parameter NLIMIT defined under keyword BASIS_SIZE. It should also be larger or equal to 
the actual size of the HO basis LDBASE, which can be larger than NLIMIT in case of degenerate 
HO states. 

PARAMETER ( NDXHRM=33 , NDYHRM=33 , NDZHRM=45 ) 

Correspond to the maximum numbers of the HO Gauss-Hermite integration nodes in three 
Cartesian directions. Should be larger or equal to the input parameters NXHERM, NYHERM, and 
NZHERM, respectively, that are defined under keyword GAUSHERMIT or determined in an optimal 
way, see keyword OPTI_GAUSS. 

PARAMETER (NDSTAT=181) 

Corresponds to the maximum numbers of the single-particle states or quasiparticle states, which 
are kept after diagonalization. Should be larger or equal to the input parameters NUMBSP(0,0) 
+ NUMBSPd.O) and NUMBSP(0,1) + NUMBSP(1,1), see keyword PHASESPACE. 

On vector machines, parameters NDBASE, NDXHRM, NDYHRM, NDZHRM, and NDSTAT should be 
odd integers in order to minimize the risk of bank memory conflicts. By the same token, 
parameter NDMAIN should be even, because it defines the array dimensions beginning with 0. 

4.2. Vectorization properties 

As discussed in Sec. 1-4.3, the code hfodd has to operate by using seven-fold nested short 
loops, and this part does not perform well in a vector processor. However, it turns out that 
the loops can be artificially made longer in such a way that the final CPU time in a vector 
processor actually becomes much shorter. All the places where this trick has been applied can 
be identified in the source file by finding the lines beginning with CVECTOR and CSCALAR, for 
example: 



CVECTOR 



DO KZ=0 , LAZOXY (NX , NY) +LAZOXY (MX , MY) 
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CSCALAR 



DO KZ=NZ+MZ,0,-2 

RESULT=RESULT+COEFOO (KZ , NZ ,MZ , 3) *T_AUXI (KZ) 
END DO 



The line beginning with CSCALAR should be made active on a scalar or superscalar machine. 
The line immediately below the line which begins with CVECTOR should be active on a vector 
machine. The results of calculation do not depend on which version of loops is activated. 



The code hfodd version (v2.40h) is written in FORTRAN- 77. However, on several platforms, 
several tested FORTRAN-77 compilers did not perform well with respect to the memory man- 
agement. In particular, they often distributed large arrays on computer's stack and heap in 
a wasteful way, and then the available memory could be insufficient for the code to execute, 
especially for a large HO basis. Therefore, crucial parts of the code have been promoted to 
FORTRAN-90, with a use of an explicit memory allocation and deallocation. In the FORTRAN 
source code provided in the file hfodd. f, all these FORTRAN-90 features are commented out 
and inactive. However, very simple modifications of the source code can easily be performed 
to transform code HFODD to FORTRAN-90. To this end, the user should run an automated 
editor script which performs the following two operations: 

1. Replace all strings "C F90" by spaces, and 

2. Remove two lines of code (or three lines of the file) after lines marked by the string: 
"C IN F90 VERSION PLEASE REMOVE THE TWO FOLLOWING LINES". 

The FORTRAN-90 version should be compiled with the "fixed form" non-standard option of 
the FORTRAN-90 compiler. 

A set of c-shell and ex-editor scripts is provided within the hfodd distribution file, and 
allows performing the above modifications of the FORTRAN source code. 

The present version of the code hfodd version (v2.40h) is the last one that works under 
FORTRAN-77; future releases will only use FORTRAN-90 programming. 

4.4. Library subroutines 
4.4.1. BLAS 

The code hfodd requires an implementation of the BLAS (Basic Linear Algebra Subprogams) 
interface for common dense vector and matrix computations. A reference implementation of 
these subroutines can be downloaded from 
|http : / / www . net lib . org/blas /bias . t gz , 

but for peak performance it is recommended that an optimized version of these routines should 
be installed. Optimized machine-specific BLAS (such as Sun Performance Library and Intel 
Math Kernel Library) are available from most hardware vendors, and there are also third-party 
optimized implementations such as GOTO and ATLAS available for various architectures. 

The BLAS subroutines are in the REAL*8/C0MPLEX*16 version, and should be compiled 
without promoting real numbers to the double precision. On the other hand, the code hfodd 
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FORTRAN-90 version 
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itself does require compilation with an option promoting to double precision. Therefore, the 
code and the BLAS package should be compiled separately, and then should be linked together. 

4.4.2. Diagonalization subroutines 

The code hfodd requires an external subroutine that diagonalizes complex hermitian matri- 
ces. Version (vl.60r) (see II) has been prepared with an interface to the NAGLIB subroutine 
f02axe, version (vl.75r) (see III) with an interface to the LAPACK subroutine zhpev, and 
version (v2.08i) (see IV) with an interface to the LAPACK subroutine zhpevx. In the present 
version (v2.40h), all these interfaces remain supported and can be activated as described in 
II-IV. However, the recommended interface is now to the LAPACK subroutine zheevr, as 
described in this section below. 

In version (v2.08i) we have implemented interface to the LAPACK subroutine zhpevx, 
which can be downloaded (with dependencies) from 

pittp: / /netlib2.cs.utk.edu/cgi-bin/netlibfiles.pl?filename=/lapack/complexl6/zhpevx.f 
This subroutine finds not all, but only the lowest eigenvectors, and hence performs calcula- 
tions in less CPU time; the gain is particularly significant for a large HO basis. Numbers of 
eigenvectors to be found are defined by the size of the phasespace, see keyword PHASESPACE. 

Subroutine zhpevx and its dependencies are in the REAL*8/C0MPLEX*16 version, and should 
be compiled without promoting real numbers to the double precision. Therefore, the code and 
the ZHPEVX package should be compiled separately, and then should be linked together. 

In order to activate the interface to the LAPACK zhpevx subroutine, the following modi- 
fications of the code hfodd (v2.40h) have to be made: 

1. Change everywhere the value of parameter I_CRAY=1 into I_CRAY=0. 

2. Change everywhere the value of parameter IZHPEV=0 into IZHPEV=2. 

3. If your compiler does not support undefined externals, or subroutines called with different 
parameters, remove calls to subroutines cgemm, f02axe, zhpev, and zheevr. 

In version (v2.40h) we have implemented interface to the LAPACK subroutine zheevr, 
which can be downloaded (with dependencies) from 

http : //www. net lib . org/ cgi-bin/netlibf iles .pl?f ilenaine=/lapack/ complexl6/zheevr . f 
This subroutine works with unpacked matrices and hence performs calculations in less CPU 
time at the expense of using a larger memory. Alternatively, the entire LAPACK library can 
be downloaded from 

http : //www . net lib . org/lapack/lapack . tgz 

Subroutine zheevr and its dependencies are in the REAL*8/C0MPLEX*16 version, and should 
be compiled without promoting real numbers to the double precision. Therefore, the code and 
the zheevr package should be compiled separately, and then should be linked together. 

In order to activate the interface to the LAPACK zheevr subroutine, the following modi- 
fications of the code hfodd (v2.40h) have to be made: 

1. Change everywhere the value of parameter I_CRAY=1 into I_CRAY=0. 

2. Change everywhere the value of parameter IZHPEV=0 into IZHPEV=3. 
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3. If your compiler and linker do not support undefined externals, or subroutines called with 
different parameters, remove calls to subroutines C'GEMM, f02axe, zhpev, and zhpevx. 

4.4.3. Matrix inversion subroutines 

The code hfodd requires external subroutines that invert complex matrices and calculate their 
determinants. In version (v2.40h) we have implemented interfaces to the LINPACK subroutines 
ZGEDI and ZGECO, which can be downloaded from 



http://www.netlib.Org/linpack/zgedi.f , and 



http://www.netlib.Org/linpack/zgeco.f , respectively, together with 



http: / / www.netlib.org/linpack/ zgefa.f 



These subroutines are in the REAL*8/C0MPLEX*16 version, and should be compiled without 
promoting real numbers to the double precision. Therefore, the code and these subroutines 
should be compiled separately, and then should be linked together. 

4.5. Printing the execution times 

After successful completion of execution, the code hfodd prints the table of CPU times spent 
in principal subroutines, as well as the total execution time. Since implementation of the 
time-checking functions widely varies between compilers and platforms, in the file hfodd. f all 
references to these functions are commented out and inactive. However, very simple modifica- 
tions of the source code can easily be performed to activate one of the following options: 

1. Replace all strings "CCPUT" by spaces in order to activate a call to the FORTRAN-90 
subroutine cpu_time, which returns CPU execution times. 

2. Replace all strings "CETIM" by spaces in order to activate a call to function etime, which 
returns CPU execution times. 

3. Replace all strings "CUNIX" by spaces in order to activate a call to the UNIX function 
TIME, which returns wall clock execution times. 

4. Replace all strings "CSAL" by spaces in order to activate a call to the SALFORD© 
FORTRAN subroutine dclock(Q), which returns wall clock execution times. 

5. Replace all strings "CF32" by spaces in order to activate a call to the MICROSOFT© 
F32 subroutine gettim, which returns wall clock execution times. 



Other time-checking functions, available to the user, can be analogously called in the subroutine 

CPUTIM. 



5. TEST RUNS 



73 



5. TEST RUNS 

Several examples of input and output files are provided with the released versions (vl.60r), 
(vl.75r), (v2.08i), and (v2.40h) of the program hfodd, or added in this user's guide, as listed 
in Table O The input data files dyl52-a.dat through dyl52-f.dat present the valid input 
items of each version and the default or recommended values of the input parameters. The 
input data files tbl51-a.dat through tbl51-d.dat are distributed with version (vl.75r), and 
provide examples of calculations with diabatic blocking. The input data files lal32-a.dat 
and lal32-b.dat are specific for version (v2.08i) and later, and illustrate various symmetry 
options of the code. The files pbl28-a.dat and snl20-a.dat are added in this user's guide 
and give examples of the spherical HF and HFB calculations, respectively. The input data 
files ge064-a.dat through ge064-c.dat and snl20-b.dat are specific for version (v2.40h) and 
later, and illustrate the AMP calculations and Lipkin-Nogami method, respectively. 



Table 6: Input data files supplied with versions (vl.60r), (vl.75r), (v2.08i), and (v2.40h) of the 
program HFODD or added in this user's guide. Each output file listed has been obtained by 
running the version of the code specific for the given input data file. 



Distribution 


Input 


Contents 


Output 


(vl.eOr) 


dyl52-a 


dat 


Only keyword EXECUTE 

Example of superdeformed HF calculations 


dyl52-a 


out 




dyl52-b.dat 


All input items of (vl.GOr), 


dyl52-b 


out 








default values of input items 






(vl.75r) 


dyl52-c 


dat 


All input items of (vl.75r), 
default values of input items 


dyl52-c 


out 




dyl52-d 


dat 


Input items added in (vl.75r), 
recommended values of input items 


dyl52-d 


out 




tbl51-a 


dat 


Example of adiabatic calculations 


tbl51-a 


out 




tbl51-b 


dat 


Example of adiabatic calculations 


tbl51-b 


out 




tbl51-c 


dat 


Example of diabatic blocking 


tbl51-c 


out 




tbl51-d 


dat 


Example of diabatic blocking 


tbl51-d 


out 


(v2.08i) 


dyl52-e 


dat 


All input items of (v2.08i), 
default values of input items 


dyl52-e 


out 




lal32-a 


dat 


Example of tilted-axis-cranking calculations 


lal32-a 


out 




lal32-b 


dat 


Example of tilted-axis-cranking calculations 


lal32-b 


out 


(v2.40h) 


dyl52-f 


dat 


All input items of (v2.40h) 
default values of input items 


dyl52-f 


out 




ge064-a 


dat 


Example of AMP calculations 


ge064-a 


out 




ge064-b . dat 


Example of AMP calculations 


ge064-b 


out 




snl20-b 


dat 


Example of Lipkin-Nogami method 


snl20-b 


out 


User's 


pbl28-a 


dat 


Example of spherical HF calculations 


pb208-a 


out 


guide 


snl20-a 


dat 


Example of spherical HFB calculations 


snl20-a 


out 
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5.1. Descriptions of the test runs 
5.1.1. Files dyl52-a.dat through dyl52-f.dat 

As explained in Section [2D) all the input data parameters of the code have their default values, 
which are assumed if the relevant input items are absent in the input data file. Each parameter 
introduced in a given version is also valid in all later versions, and its default value remains 
unchanged. Later versions introduce new parameters, but their default values are such that the 
new features are not active. Thus, the results of running a later version are identical to those 
of running an earlier version, provided that the values of those new parameters are left at the 
default. In particular, running each version with its default data set corresponds to identical 
conditions and gives the same solution. All the input items valid for a given version are also 
recognized and read by the subsequent ones. In general, running a later version with an input 
data file prepared for an earlier version leads to the same solution. However, the structure of 
the output file may be modified from version to version; in particular because later versions 
may calculate more observables. 

Running the code with the default values results in performing 50 iterations for the super- 
deformed state in ^^^Dy at the angular frequency of hu = 0.5 MeV. The pairing correlations 
are not taken into account. The HO basis is defined by an axially symmetric (with respect 
to the z axis) nuclear surface with physical deformation typical for superdeformed bands in 
the A ~ 150 region. The parity, y-signature, x-simplex^ (and their products) are imposed as 
conserved symmetries, while the time-reversal symmetry is broken in order to describe rota- 
tion (about the y axis). The run starts from the Nilsson potential, and an appropriate weak 
constraint on the quadrupole moment (Q20) is used to ensure that the iterations fall into the 
superdeformed minimum. 

The input data files dyl52-a.dat through dyl52-f.dat define either the default or rec- 
ommended values of the input parameters for the released versions of the code. They contain 
either all or some of the valid input items, as described in Table El File dyl52-a.dat contains 
only keyword EXECUTE, and when used as input for any version of the code it results in per- 
forming calculations with the default data set. Files dyl52-b.dat, dyl52-c.dat, dyl52e-dat, 
and dyl52f-dat contain all the valid input items of versions (vl.GOr), (vl.75r), (v2.08i), and 
(v2.40h), respectively, and the input data are identical to the default values. File dyl52-d.dat 
contains only the input items introduced in version (vl.75r), and the corresponding new param- 
eters are set to recommended values, which differ from the default ones in the following way. 
The default values are such that the code works identically to version (vl.60r). The parameters 
EPSITE, NULAST, EPSPNG, and NUCHAO, control the termination of the iterations. By default, 
they are set so that always the number of iterations specified by NO ITER is performed. The 
recommended values define an optimal level, at which convergence, "chaotic" divergence, or 
"ping-pong" divergence are detected to stop the iterations earlier. The differences between the 
default and recommended values of the variables governing the diabatic blocking are irrelevant 
because the diabatic blocking is disabled by default. Also by default, the Coulomb file is neither 
read from nor written to disc. With the recommended values of the parameters ICOULI and 
ICOULO, its reading and writing is handled automatically, which saves the CPU time. 

Files dyl52-a.dat through dyl52-f .dat serve only as suitable patterns to modify the input 
parameters. However, it should be a good practice to include in the input data file only those 
input items which modify the input parameters with respect to the default values. In this 
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J I I I I I I II I I I I I I L 



0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 
Rotational frequency (MeV) 

Figure 5: Negative-parity neutron single-particle Routliians in ^^^Tb calculated for the 
[761]3/2(r=+i) (a) and [514]9/2(r=+i) (b) neutron diabatic configuration. Solid and dashed 
curves denote the r=+i and r=-i signatures, respectively. The arrows denote the angular 
frequencies where the converged solutions near the crossing points could not be found. 

way, the danger of accidentally modifying some intricate numerical parameters of the code is 
minimized. 

5.1.2. Files tbl51-a.dat through tbl51-d.dat 

Files tbl51-a.dat through tbl51-d.dat provide data for calculations in which single-particle 
configurations corresponding to two crossing bands in ^^^Tb [12] are followed either adiabatically 
or diabatically. The pairing correlations are not taken into account. The parity, y-signature, 
x-simplex^ (and their products) are imposed as conserved symmetries, while the time-reversal 
symmetry is broken in order to describe rotation (about the y axis). The angular frequency 
range from hu = 0.50 to 0.80 MeV is scanned with the step of 0.05 MeV. Calculations at 
hu = 0.50 are started from the Nilsson potential, and iterations for each subsequent frequency 
are restarted from the solution obtained at the preceding frequency. The two considered bands 
are based on excited physical configurations, in which the last neutron occupies either the 
[514]9/2(r=+i) or [761]3/2(r=+i) Nilsson level. At hoo = 0.50 MeV the level [514]9/2(r=+i) 
is lower in energy, and has number 23 --|- (23rd state of 7r=-l and r=+i), while the level 
[761]3/2(r=-|-i) has number 24--I-. These two levels cross in the concerned frequency interval 
so that at hu = 0.80 MeV the situation is the opposite. This is illustrated in Figure [5], which 
is explained in detail below. 
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In the files tbl51-a.dat and tbl51-b.dat it is required that always the level 23-+ or 24- 
+ be occupied, respectively, which amounts to following the configuration adiabatically. If, say, 
the level 23-+ is occupied, the solutions at low and high frequencies correspond to physical 
configurations [514]9/2(r=+z) and [761]3/2(r=+z), respectively, while in the vicinity of the 
crossing no solution can be found, because the so-called "ping-pong" divergence occurs. This 
means that the physical states [514]9/2(r=+z) and [761]3/2(r=+«) change their energy order 
every second iteration, so that each of them once becomes the level 23-+. Thus, putting the 
neutron always on the level 23-+ results in occupying alternately the very different physical 
states [514]9/2(r=+i) and [761]3/2(r=+2), which cannot lead to convergence. 

The code detects the "ping-pong" divergence, stops the iterations, and prints a table that 
helps identify the crossing orbitals. Such a table is reproduced in Table[7]for the 23- + adiabatic 
configuration at hu = 0.65 MeV. For both protons and neutrons, in each parity-signature block, 
and separately for the particle and hole states, the code finds the level for which the difference 
in the angular momentum or intrinsic spin alignment (on the y axis) between the last and last 
but one iteration is the largest. The indices of those levels in the parity-signature blocks, and 
the absolute values of the differences are printed in the upper part of the table, independently 
for the angular momentum and intrinsic spin. If in some symmetry block, the indices of the 
particle and hole states found in this way differ by one, then it is probable that these two 
levels cross and cause the "ping-pong" divergence. When such pairs exist, they are listed in the 
lower part of the table. In the given example, the pair of states 23-+ and 24-+ is properly 
recognized on the basis of both the angular momentum and intrinsic spin. The crossing orbitals, 
[514]9/2(r=+i) and [761]3/2(r=+z), have markedly different alignments of about —0.46^ and 
+2.23, respectively, which corresponds to the values 2.820 and 2.807/i given in the table. 

If one wishes to follow, say, the configuration [761]3/2(r=+i) diabatically, one should request 
that, of the levels 23-+ and 24-+, the one with physical properties of the state [761]3/2(r=+z) 
be occupied. One can chose, for instance, the state with higher angular momentum alignment, 
independent of whether it is lower or higher in energy. In order to diabatically follow the 
configuration [514]9/2(r=+z), one can chose the state with lower alignment. The corresponding 
runs are defined in the files tbl51-c.dat and tb-151-d.dat. Indeed, in each case solutions 
corresponding to the same physical configuration is found at each angular frequency. Moreover, 
the "ping-pong" divergence does not occur, even very close to the crossing frequency. Figure [5] 
shows the negative-parity neutron single-particle Routhians obtained from the calculations of 
files tbl51-c.dat and tb-151-d.dat. 

5.1.3. Files lal32-a.dat and lal32-b.dat 

File lal32-a.dat defines two consecutive runs, which lead to a planar solution [13] in the 
triaxially deformed nucleus ^'^^La for hu = 0.25 MeV. The pairing correlations are not taken 
into account. Since the calculations involve the tilted-axis cranking, a spherically invariant 
HO basis (with the three HO frequencies equal and comprising only entire HO shells) is used 
in order not to favor any orientation of the solution in space. The first run is performed at 
hu = with the parity, ^/-signature, x-simplex^ (and their products) conserved, and with the 
time-reversal broken in order to describe odd numbers of neutrons and protons. Due to the 
presence of the odd nucleons, even for huj=0 one obtains a non-zero angular momentum, which 
is oriented along the y axis. The run starts from the Nilsson potential, and in 94 iterations finds 
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Table 7: Example of the output printed when the "ping-pong" divergence is found, see text. 

* * 

* TABLE BELOW GIVES THE MAXIMUM ABSOLUTE VALUES OF DIFFERENCES OF ALIGNMENT * 

* BETWEEN THE LAST TWO ITERATIONS. IT MAY SERVE AS A GUIDE TO SELECT STATES * 

* WHICH EXCHANGE WAVE FUNCTIONS IN EVERY SECOND ITERATION. SUCH STATES MAY * 

* BE CANDIDATES FOR A CALCULATION WITH DIABATICALLY BLOCKED CONFIGURATIONS * 

* * 

* * * 

* BLOCKS * ANGULAR MOMENTUM ALIGNMENT INTRINSIC SPIN ALIGNMENT * 

* * * 

* * PARTICLES HOLES PARTICLES HOLES * 

* ISO. PAR. SIG. * INDEX VALUE INDEX VALUE INDEX VALUE INDEX VALUE * 

* * * 



* 








* 


























* 


* 


NEUT 


+ 


+ 


* 


38 


0. 


.221 


13 


0. 


.018 


24 


0. 


.023 


13 


0. 


.016 


* 


* 


NEUT 


+ 




* 


28 


0. 


.163 


7 


0. 


.041 


33 


0. 


.014 


7 


0. 


.032 


* 


* 


NEUT 






* 


39 


0. 


.095 


13 


0. 


.198 


44 


0, 


.015 


14 


0. 


.150 


* 


* 


NEUT 




+ 


* 


24 


2, 


.820 


23 


2. 


.807 


24 


0. 


.164 


23 


0. 


.167 


* 


* 


PROT 


+ 


+ 


* 


24 


0. 


.255 


13 


0. 


.149 


25 


0. 


.129 


14 


0. 


.122 


* 


* 


PROT 


+ 




* 


27 


1. 


.733 


13 


0. 


.062 


21 


0. 


.319 


8 


0. 


.020 


* 


* 


PROT 






* 


32 


0. 


.119 


10 


0. 


.530 


32 


0. 


.029 


10 


0. 


.065 


* 


* 


PROT 




+ 


* 


24 


2. 


.649 


15 


0. 


.082 


24 


0. 


.153 


10 


0. 


.030 


* 


* 








* 


























* 



* * 

* FOLLOWS THE LIST OF CANDIDATE CONFIGURATIONS FOR THE DIABATIC BLOCKING * 

* * 

* NEUT - + * 24 2.820 23 2.807 * 

* NEUT - + * 24 0.164 23 0.167 * 

a triaxial solution, whose shape is forced by using appropriate constraints for the quadrupole 
moments {Q20) and {Q22)- The constraints are set so that this solution nearly corresponds to 
the unconstrained minimum, and are used only to prevent the iterations from falling into other 
local minima. These constraints are chosen so that the short, medium, and long axes of the 
mass distribution coincide with the y, x, and z axes, respectively. For the concerned solution, 
the single-particle configuration defined in the input data file in terms of numbers of particles 
in the parity-signature blocks physically corresponds to the odd proton particle occupying the 
lowest substate of the /iii/2 orbital and the odd neutron hole left on the highest /iii/2 substate. 

The second run starts from the solution obtained in the first run, and performs calculations 
with all symmetries broken except parity, with multipole constraints released, and for the 
angular frequency of hu = 0.25 MeV. The single-particle configuration specified in terms of 
parity blocks is obtained from the configuration for the first run simply by adding the numbers 
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of particles in the two signature sub-blocks of each parity block. All the three components of 
the initial angular frequency vector are non-zero, hux = fiUy = huiz = (0.25/\/3) MeV, which 
breaks all symmetries except parity in the first iteration, thus making the calculations otherwise 
symmetry-unrestricted. If a symmetry was not broken in the first iteration, it would artificially 
remain conserved as governed by the theorem about self-consistent symmetries 

The Kerman-Onishi condition [T3] states that the angular momentum and angular frequency 
vectors must be parallel in self-consistent solutions. Standard iteration, with IMOVAX=0, re- 
quires several thousand iterations for the nucleus to turn appropriately in space to achieve 
the alignment of these vectors below the angle of 0.01°. On the other hand, iteration with 
IM0VAX=1, as shown in the example of file lal32-a.dat, achieves the alignment of 0.000001° 
within 92 iterations. In the solution resulting from the second run, the angular momentum 
vector has non-zero components only on the y and z axes of the intrinsic frame, which are the 
short and long axes of the triaxial mass distribution. The calculations of the file lal32-a.dat 
are performed in two steps to get a better control over what is going on and to avoid divergence 
or falling into undesired minima. 

The input data provided in the file lal32-b . dat are similar to those of the file lal32-a. dat, 
except that the second run is performed with the x-simplex"^ conserved. In this way, the 
angular momentum and angular frequency vectors are confined within the plane y-z, which 
is also a principal plane of the mass distribution. After the first run, this plane contains 
the short and long axes. The components of the initial angular frequency vector are set to 
hujy = hu, = (0.25/v/2)MeV and hux=0. Since the planar solution in question corresponds 
to non-equal projections of the angular momentum on the short and long axes, in order to 
satisfy the Kerman-Onishi condition, the nucleus still has to turn in space about the x axis. 
For IMOVAX=0, this process is also very slow, while for IM0VAX=1, as in file lal32-b.dat, 
it takes again only 92 iterations. The x-simplex-^-imposing solution obtained by using the 
file lal32-b.dat differs from the symmetry-unrestricted result of the file lal32-a.dat only by 
orientation in space; apart from that, the values of all observables are identical. This shows that 
in the considered self-consistent solution the x-simplex-^ symmetry is spontaneously conserved. 

5.1.4. File pb208-a.dat 

The input data file pb208-a.dat defines two consecutive runs to perform precise calculations 
for the spherical ground state of the doubly-magic nucleus ^osp]-, large basis of 26 HO 
shells. The pairing correlations are not taken into account. Conservation of all the symmetries 
considered in the code (the group -Dj^) is imposed. For the spherical shape, the specified 
parity-signature configuration corresponds to occupation of the entire magic shells, up to the 
proton and neutron magic numbers 82 and 126, respectively. A spherical basis is used, i.e., 
with all the three characteristic HO frequencies equal {huJx=fiujy=huJz=fiu!o) and comprising 
only the entire HO shells. The so-called physical value of hcuQ = / x 41 MeV/A^^^ with / = 1.2 
is taken, following the default settings of the code. To speed up the calculations, the program 
HFODD uses a diagonalization algorithm that can reduce the computational effort by finding 
less eigenstates than the matrix size (the lowest ones). In the present example, in each charge- 
simplex block only one extra state is calculated as compared to the number of occupied states in 
that block, as governed by the input parameter NUMBSP. This is far less than the code's default, 
and such a way of proceeding is possible only because no particle-hole excitations into higher 
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states are required. Also for the sake of saving time, a solution in a smaU basis of A'"o=12 HO 
shells is obtained first, starting from the Nilsson potential, and then serves as a good starting 
point for the time-consuming calculations with Nq=26. The number of states in the spherical 
basis comprising the HO shells from to No equals M={Nq + l)(iVo + 2)(A^o + 3)/6, and hence 
the numbers of the HO basis states (for each projection of the intrinsic spin) equal M=455 and 
M=3654, respectively for iVo=12 and Nq=2Q. Since the code hfodd can restart iterations only 
if the number of Gauss-Hermite integration nodes is preserved, both runs are performed with 
54 points in each of the Cartesian directions, which is excessive for A''o=12, but optimal for 
Nq=2Q. For Nq=2Q, the 54 points ensure exact integration of all matrix elements apart from 
those of the Coulomb interaction and density-dependent term. 

The energies of ^osp^ obtained from calculations with 7Vo=12 and A'"o=26 equal —1630.638 
and —1635.531 MeV, respectively. The latter value is very close to the result of —1635.666 MeV 
obtained by using a very precise one-dimensional spherical code working in spatial coordinates 
jl5] . Although the Nq=12 HO basis is evidently too small for a precise determination of the 
energy, other observables may not be so sensitive to the HO basis size. For instance, the mass 
root-mean-square radii obtained for A''o=12 and A^o=26 equal 5.5540 and 5.5554 fm, respectively, 
the result of coordinate-space calculations being 5.5550 fm. 

5.1.5. File snl20-a.dat 

The input data file snl20-a.dat provides an example of the HFB calculations for the spherical 
ground state in ^^°Sn. All symmetries considered in the code (the group -Dj/i) are imposed, 
except the particle-number symmetry, which is broken in order to describe the pairing corre- 
lations. A spherical basis is used. The pairing correlations are included within the full HFB 
method, and the zero-range density-dependent interaction given by the pairing form-factor 
([TT]) is used in the particle-particle channel. The particle and pairing densities are obtained by 
summing up contributions from the quasiparticle states with the equivalent-spectrum energies 
not exceeding emax=60 MeV, which is the code's default pairing window. In the pairing form- 
factor (fTTj) . the power a=l and the saturation density po=0.32fm~^ are taken. The interaction 
strength, Vb=-285.88 MeVfm^, was adjusted for the given cutoff energy, emax? to reproduce 
the value of A7v=1.245MeV for the average neutron pairing gap in ^^°Sn. The same values of 
a, po and Vq are taken for both neutrons and protons. The Fermi energy of —8 MeV and the 
pairing gap of 1 MeV are used in the first iteration to construct the HFB Hamiltonian, while 
the initial mean field is given by the Nilsson potential. 

Since ^^°Sn is a proton-magic nucleus, the proton pairing decreases during the iteration and 
vanishes in the HFB solution. As expected, one obtains a significant thickness of the neutron 
skin (difference between the neutron and proton root-mean-square radii), which in the present 
case equals 0.14fm. 

5.1.6. Files ge064-a.dat and ge064-b.dat 

File ge064-a.dat contains the input data that allow for determining the rotating triaxial 
state in the nucleus ^^Ge and then performing a 3D AMP of this state. This input data file 
comprises four runs of the code: (i) starting from the Nilsson initial potential, the code performs 
20 iterations with constrained values of quadrupole moments Q20 and Q22, such that the initial 
deformation and orientation of the nucleus is determined, (ii) after releasing the constraints, the 
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triaxial minimum is converged, (iii) the rotating solution for / = 6 is converged by setting the 
angular frequency of huj = 0.575 MeV, and (iv) the AMP of the rotating solution is performed. 
Small numbers of the Gauss-Tchebyschev (10) and Gauss- Legendre nodes (10), see Section VI- 
2.1, which are used in this run, do not allow for a precise AMP resolution of high angular 
momenta. This input data file is only meant to provide an example of a rapid calculation. 

File ge064-b . dat contains the input data that allow for a correct AMP, with higher numbers 
of the Gauss-Tchebyschev (40) and Gauss-Legendre nodes (40), but its execution requires a 
CPU time which is 4^=64 times longer (about 200h). Alternatively, one can run in parallel 
40 jobs by executing the input data files ge064-c.dat with characters NUASTA replaced by 
integers from 1 to 40. 

5.1.7. File snl20-b.dat 

File snl20-b.dat contains the input data that allow for determining the ground state in the 
nucleus ^^°Sn with the Lipkin-Nogami corrections taken into account. 
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NYHERM, 


46, 
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BLOCKSIM_N, 37 
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EVE_ADD_TS, 18 
EVE_SCA_PM, 17 
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ODD_SCA_PM, 18 
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PHPARI_NEU, 38 
PHPARI_PRO, 38 
PHSIGN_NEU, 41 
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PHSIMP_NEU, 36 
PHSIMP_PRO, 36 
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PRINT-INTR, 55 
PRINT-ITER, 55 
PRINT-MOME, 55 
PRINT_AMP, 56 
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REPLAYFILE, 58 
REPLIPFILE, 61 
REPYUKFILE, 60 
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SLOWLIPKIN, 33 
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SPIN_ORBIT, 14, 15 
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TSIMPLEX3D, 27 
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WOODSAFILE, 59 
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CPUTIM, 72 

DCLOCK, 72 
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F02AXE, 71, 72 
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PARAMS, 14 
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